{ "cells": [ { "cell_type": "code", "execution_count": 55, "id": "optimum-atlanta", "metadata": {}, "outputs": [], "source": [ "# In this Jupyter notebook we will fit a nonlinear function to a set of experimental\n", "# data. The fit will be weighted by the measurement uncertainties.\n", "\n", "# In this example, we take data from a series LRC circuit. For the y-axis\n", "# we will have the ratio of the magntidue of the voltage across the\n", "# resistance to that supplied by the function generator driving the\n", "# circuit. On the x-axis we will have frequency.\n", "\n", "# The process will be very similar to the process used to fit a linear function to\n", "# a dataset. The main difference will be in the definition of the model function\n", "# that we are fitting to.\n", "\n", "# Import the NumPy, Matplotlib and SciPy modules\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.optimize import curve_fit" ] }, { "cell_type": "code", "execution_count": 56, "id": "chronic-pricing", "metadata": {}, "outputs": [], "source": [ "# Enter the data as arrays.\n", "Vratio= np.array([0.198e-1, 0.67e-1, .117, .185, .331, .450, .573, .689, .718,\\\n", " .714, .704, .670, .631, .549, .412, .318, .249, .186, .138])\n", "errVratio= np.array([0.1e-2, 0.2e-2, 0.3e-2, 0.3e-2, 0.5e-2, 0.6e-2, 0.7e-2,\\\n", " 0.8e-2, 0.8e-2, 0.8e-2, 0.8e-2, 0.8e-2, 0.8e-2, 0.7e-2,\\\n", " 0.6e-2, 0.4e-2, 0.5e-2, 0.3e-2, 0.3e-2])\n", "f = np.array([3000, 10000, 15000, 20000, 25000, 28000, 30000, 32000, 33000,\\\n", " 34000, 35000, 36000, 37000, 39000, 43000, 47000, 52000, 60000,\\\n", " 70000])" ] }, { "cell_type": "code", "execution_count": 57, "id": "talented-trailer", "metadata": {}, "outputs": [], "source": [ "# Calculate angular frequency.\n", "omega = (2*np.pi*f)" ] }, { "cell_type": "code", "execution_count": 58, "id": "municipal-chinese", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAhZ0lEQVR4nO3de3xcdZ3/8debSbGCWy5SFdpAC8vFwk/SEltRl4XKpYXUyIIPuVSrrttfFVZdcRV+Ig3t+sPL6k9ccWth2dUtFy9cLJG2SBVdL902kICltFi5NUVpuQgIujTh8/vjnNBJOplMmjmZTOb9fDzmkTnfc+bMZ75N85lzvud8vooIzMysdu1R6QDMzKyynAjMzGqcE4GZWY1zIjAzq3FOBGZmNa6u0gEM1gEHHBCTJk2qdBhmZlXl7rvvfjIixhdaV3WJYNKkSbS1tVU6DDOzqiLp0f7W+dSQmVmNcyIwM6txTgRmZjXOicDMrMY5EZiZ1TgnAjOzGudEYGZW46ruPgKzLHR3w4oV0N4OU6fC7NmQy1U6KrPh4SMCq1ktLS1IQsoxdtxK5l20lRXP38u8i7YydtxKpBySaGlpqXSoZplStU1M09jYGL6z2MqptRUuWPgsc9csJzcm6N4h/nXKTJ7ZfA3QBbQDK4CXAVi4cKGTg1UdSXdHRGOhdT41ZDWvvR3qT91CbkzypUh7QG7M3rz2iI9wZPPv6Fz9SaZM2JvWW+p8ushGJScCq3lTp8I1C+vpXtRObkzw4G0TIMT/Xr8yPUK4l2UzmmltHUcu53EEG30yHSOQNEvSJkmbJV1cYP0/SupIH+sldUvaP8uYzPqaPRumTNibbx5zMqs/1cAdH3krR7xz5xFCbkxw0MyHOfu8TR5HsFEpszECSTngQeAUoBNYB5wbERv62X4O8A8RMbPYfj1GYOXS0tLC5Zdfni7tAcwGGoAx7H/Eh1iw/sevjBlcfdRZ7PmaHcxr2zmOsGxGM1ctGkdTU8U+glnJio0RZHlEMB3YHBEPRcRLwI1Ac5HtzwVuyDAes15aWlqIiPTRzcKFjcDngEU81/lrlhw9kzs/dSxLjp7JU1u2cvBpj/Q6Sqg/7TE6Oir5CczKI8tEMAHYkrfcmbbtQtJewCzgpgzjMStqZ2Lo5s/PzeJbX5nA6eOO5VtfmcDN3zucrasn0b1DAHTvEFtWHUxDQ2VjNiuHLBOBCrT1dx5qDvCLiHi64I6k+ZLaJLVt3769bAFa7dh5z4DS8/pNSJ9Nf+ZeWddzrj+Xg6YmuPTS5GdTU+9xhGUzmpkycS9mz67s5zIrhyzHCI4HWiLitHT5EoCIuKLAtrcA34uI6wfar8cIbCi6u6HpzC7WPbiNI9/5JJ2rDyl6aWj/4wgd+N4CqybFxgiyTAR1JIPF7wC2kgwWnxcR9/fZbh/gYaA+Il4YaL9OBDYUhW4e86Cv1YKKDBZHRBdwIbAKeAD4bkTcL2mBpAV5m54J3FFKEjAbqr43j3nQ1yzj+wgi4vaIOCIiDouIz6VtSyJiSd42/xER52QZh1mPqVNhyx31HvQ1y+Oic1YTegaL58zJ8fjGX/W6NPTxTb9kzpzB3xjW3Z2calq8OPnZ3Z1d/GZZctE5qzk9Jac7OqChYXClInYOHu9B3V4/ZNzE/8WRzU+y6QcH8Fznr+l68QzgZQ8e24hTkcHirDgR2EjgQWerNpW6s9hs1PKgs40mTgRmu8GDzjaaOBGYDUIWg85mleYxArPdNJRBZ7Ph5hnKzDLQU4/Ig8NW7XxqyMysxjkRmJnVOCcCM7Ma5zECszLqGUD2BPdWTXxEYDZEOye9yTF23EpPcG9Vx5ePmpWJy07YSOYSE2bDwGUnrFo5EZiVictOWLVyIjAbIpedsGrnMQKzMnLZCRupXGLCbJi47IRVo0xPDUmaJWmTpM2SLu5nmxMldUi6X9JPs4zHzMx2ldkRgaQccBVwCtAJrJO0PCI25G2zL/ANYFZEPCbpdVnFY2ZmhWV5RDAd2BwRD0XES8CNQHOfbc4Dbo6IxwAiYluG8ZiZWQFZJoIJwJa85c60Ld8RwH6S7pJ0t6T3FdqRpPmS2iS1bd++PaNwzcxqU5aJQAXa+l6iVAccB5wBnAZ8VtIRu7woYmlENEZE4/jx48sfqZlZDcvyqqFOoD5veSLweIFtnoyIF4AXJP0MOBZ4MMO4zMwsT5ZHBOuAwyVNlrQncA6wvM82PwD+SlKdpL2AGcADGcZkZmZ9ZHZEEBFdki4EVgE54NqIuF/SgnT9koh4QNJK4D7gZeCaiFifVUxmZrYr31lso4rnAzArzNVHbVTzfABmQ+MjAhs1PB+AWf98RGA1wfMBmO0eJwIbNTwfgNnucSKwquf5AMyGxmMENqp4PgCzwjwfgdWMapgPwJe42kjjU0Nmw8CXuNpI5lNDZsPIl7hapfjyUbMRwpe42kjkRGA2jHyJq41ETgRmw8CXuNpI5jECs2HmS1ytEnz5qNkIUg2XuFpt8akhM7Ma50RgZlbjnAjMzGqcE4GZWY1zIjAzq3GZJgJJsyRtkrRZ0sUF1p8o6VlJHenjsizjMTOzXWV2+aikHHAVcArQCayTtDwiNvTZ9L8iwhfSmZlVSJZHBNOBzRHxUES8BNwINGf4fmZmthuyTAQTgC15y51pW1/HS7pX0gpJRxfakaT5ktoktW3fvj2LWM3MalaWiUAF2vrWs7gHOCQijgX+Bbi10I4iYmlENEZE4/jx48sbpZlZjcsyEXQC9XnLE4HH8zeIiOci4o/p89uBMZIOyDAmMzPrI8tEsA44XNJkSXsC5wDL8zeQ9AZJSp9PT+N5KsOYzMysj8yuGoqILkkXAquAHHBtRNwvaUG6fglwNvBhSV3An4BzotrKoZqZVbmSylBLeidwQrr404i4LdOoinAZajOzwRvSVJWSrgA+BmxIHx9N28zMbBQo5dTQGUBDRLwMIOlbQDtwSZaBmZnZ8Ch1sHjfvOf7ZBCHmZlVSClHBFcA7ZJ+QnJvwAn4aMDMbNQYMBFExA2S7gLeTJIIPh0Rv886MDMzGx79JgJJR0XERknT0qbO9OdBkg6KiHuyD89sp55J39vbYepUT/puVi7Fjgg+AcwHvlxgXQAzM4nIrIDubmg6s4sNW1+g/tQtXLOwnilL96b1ljonA7Mh6newOCLmp09nR8RJ+Q/g9OEJz2pdS0sLkqira2LtpieYu2Y5M6+4h7lrlrN24xPU1TUhiZaWlkqHala1Srlq6JcltpmVXUtLCxHBokWtvPFvniE3JrkBMjcmeONZz7B4cSsR4URgNgT9JoK0DtBxwKslTZU0LX2cCOw1XAGaQTImsOWOerp3JEVtu3eILasOpqGhsnGZjQbFjghOA/6ZpGroV0jGCr5MMnbwf7IPzWznqaE5c3I8vvFXLDl6Jnd+6liWHD2Txzf9kjlzcj41ZDZEA9YaknRWRNw0TPEMyLWGalfPVUMdHdDQ4KuGzAajWK2hUu4juEnSGcDRwNi89kXlC9FsYLkcNDUlDzMrn1KKzi0B3gP8PckNZe8GDsk4LjMzGyalXDX01oh4H/BMRFwOHE/vmcfMzKyKlZII/pT+fFHSQcAOYHJ2IZmZ2XAqpehcq6R9gS+RTDYfwDVZBmVmA3PJDSuXAY8IImJxRPwhvXLoEOAo4POZR2Zmu+i5nFbKMXbcSuZdtJUVz9/LvIu2MnbcSiRfTmuDV/TyUUkTgAOB+yLiJUmvAz4OvD8iDhqeEHvz5aNm0NoKFyx8lrlrlpMbE3TvEMtmNHPVonG+qsoK2q2pKiV9HOgA/gVYI2ke8ADwauC4Et94lqRNkjZLurjIdm+W1C3p7FL2a1br2tuh/tQtvUpu1J/2GB0dlY3LqlOxU0PzgSMj4njgXcDVwBkR8Q8R8buBdiwpB1wFzAamAOdKmtLPdl8AVg0+fLPa5JIbVk7FEsGfI+JpgIh4DHgwItYMYt/Tgc0R8VBEvATcCDQX2O7vgZuAbYPYt1lNcskNy0K/YwSStpH88e5xTv5yRHy06I6T0zyzIuJD6fJ7gRkRcWHeNhOA60nmNvg3oDUivl9gX/NJjlA4+OCDj3v00UdL+nBmo5lLbthg7G6JiX/ss3z3YN+3QFvfrPNVkqkvu6VCm6cvilgKLIVksHiQcZiNSi65YeXSbyKIiG8Ncd+d9L4DeSLweJ9tGoEb0yRwAHC6pK6IuHWI721mZiUq5Yay3bUOOFzSZGAryaml8/I3iIhX7lCW9B8kp4ZuzTAmMzPrI7NEEBFdki4kuRooB1wbEfdLWpCuX5LVe5uZWemyPCIgIm4Hbu/TVjABRMT7s4zFzMwKK6UM9RGSVktany6/SdKl2YdmZmbDoZTqo1cDl5BUHSUi7iM5329mZqNAKYlgr4hY26etK4tgzMxs+JWSCJ6UdBjpPQDpjWIDlpgwM7PqUMpg8QUkN3MdJWkr8DAwN9OozMxs2JQyef1DwMmS9gb2iIjnsw/LzMyGy4CJQNIn+iwDPAvcHREd2YRlZmbDpZQxgkZgATAhfcwHTgSulvSp7EIzM7PhUMoYwWuBaRHxRwBJC4HvAyeQFKL7YnbhmZlZ1ko5IjgYeClveQdwSET8CfifTKIyM7NhU8oRwfUkU1X+IF2eA9yQDh5vyCwyMzMbFqVcNbRY0grgbSRzDCyIiJ7Z48/PMjgbnXomVGlvT6Zc9IQqZpVVyqkh0j/8NwA3A9skHZxpVDbq9EyxKOUYO24l8y7ayorn72XeRVsZO24lkqdYNKuUfqeqfGUD6Z3Al4GDSOYVPhjYGBFHZx/erhobG6OtrW3gDW1Eam2FCxY+y9w1y8mNCbp3iGUzmrlq0TjPtGWWoWJTVZZyRLAYeAvJ5PWTgZOBX5QxPqsh7e1Qf+oWcmOSLyC5MUH9aY/R0VHZuGpZd3eSoBcvTn52d1c6IhtupSSCHRHxFLCHpD0i4idAQ7Zh2WjTc2rossuaeODm/ejekcxR3b1DPHDTfnz2s00+NVQB3d3QdGYXFyx8llUvrueChc/SdGaXk0GNKSUR/EHSa4CfAddJuhJXH7VBamlpISLo6mpl+pGvZ9mMZn58yTSWzWhm+lGvp6urlYhwIhgmPYm5rq6JtZueYO6a5cy84h7mrlnO2o1PUFfnxFxLShkj2Bv4E0nSOB/YB1gWEU9nH96uPEZQ/XquGurogIYGXzVUSYsXw6oX1zPzinteafvxJdOYtfcxXOrpp0aVoY4RXBYRL0dEV0R8KyK+Bny6vCFaLcnloKkJLr00+ekkMPx8qs7ylXJEcE9ETOvTdl9EvGnAnUuzgCtJJq+/JiI+32d9M8lg9Mskp5s+HhE/L7ZPHxGYlU/PGMGGzhepP+0xtqw6mCkT96L1ljon6FGm2BFBvzeUSfow8BHgUEn35a36C0q4akhSDrgKOAXoBNZJWh4R+XcjrwaWR0RIehPwXeCogfZtZuWRy0HrLXWsWDGOjo5jaFjkU3W1qNidxdcDK4ArgIvz2p8vcXxgOrA5nc8ASTcCzeSVpegpZJfam3QWNDMbPj2n6nwfR+0qNkaQA54jmaHs+bwHkvYvYd8TgC15y51pWy+SzpS0Efgh8MFCO5I0X1KbpLbt27eX8NZmZlaqYkcEd7PzG7r6rAvg0AH23fc1Pa/r3RBxC3CLpBNIxgtOLrDNUpLpMmlsbPRRg5lZGfWbCNK7iIeiE6jPW54IPF7k/X4m6TBJB0TEk0N8bzMzK1EpZah76g2dkC7eFRGtJbxsHXC4pMnAVuAc4Lw++/1L4LfpYPE0YE/gqVKDNzOzoStlzuLPA28GrkubPibpbRFxSbHXRUSXpAuBVSTjDddGxP2SFqTrlwBnAe+TtIPkprX3xEDXs5qZWVmVch/BfUBDRLycLueA9lLuI8iC7yMwMxu8od5ZDLBv3vN9hhyRmZmNGKWMEVwBtEv6CcmVQCcARU8LmZlZ9Sh2Z/HXgesj4gZJd5GMEwj4dET8fpjiMzOzjBU7IvgN8GVJBwLfAW6IiI5hicrMzIZNv2MEEXFlRBwP/DXwNPDvkh6QdJmkI4YtQjMzy9SAg8UR8WhEfCEippLcB3Am8EDmkZmZ2bAo5T6CMcAskhvC3gH8FLg847jMbBTomYSovR2mTnVl05Gq3yMCSadIupakVMR84HbgsIh4T0TcOkzxmVmV6Zn0RsoxdtxK5l20lRXP38u8i7YydtxKpJwnvRlh+r2hLL1c9HrgpkpNS1mIbygzqw6trXDBwmeZu2Y5uTFB9w6xbEYzVy0a55LXFbBbN5RFxEkRcfVISgJmVj3a26H+1C3kxiRfNnNjgvrTHqOjo7Jx2a5KvbPYzKwkng+5+gxYa2ik8akhs+rg+ZBHlt2as9jMbCg8H3L1cCIws8x4PuTq4DECM7Ma50RgZlbjnAjMzGqcE4GZWY1zIjAzq3GZJgJJsyRtkrRZ0sUF1p8v6b708UtJx2YZj5mZ7SqzRJBOcn8VMBuYApwraUqfzR4G/joi3gQsBpZmFY+ZmRWW5RHBdGBzRDwUES8BNwLN+RtExC8j4pl0cQ0wMcN4zMysgCwTwQRgS95yZ9rWn78FVhRaIWm+pDZJbdu3by9jiGZmluWdxSrQVrCwkaSTSBLB2wutj4ilpKeNGhsbq6s4UhXzpCI22vh3urAsE0EnUJ+3PBF4vO9Gkt4EXAPMjoinMozHBuGVgmFbX6D+1C1cs7CeKUv3dsEwq1r+ne5flqeG1gGHS5osaU+SqS6X528g6WDgZuC9EfFghrHYIK1YARu2vsDcNcuZecU9zF2znA2dL7Ki4Mk7s5Grpyx2XV0Tazc90et3eu3GJ6irc1nszBJBRHQBFwKrSCa7/25E3C9pgaQF6WaXAa8FviGpQ5LrS1dYz3+aOXM+y4R3PNJrUpEJJz/MnDmX1vx/GqsuLS0tRASLFrXyxr95ptfv9BvPeobFi1uJiJr+nc70PoKIuD0ijoiIwyLic2nbkohYkj7/UETsFxEN6aNgrWyrhHY2/eCAXpOKbLr1AKCjolGZDZYnyhmYJ6axgjypiI02tf47XWxiGicC61fPFRYdHdDQ4CssrPrV8u+0E4GZWY0rlghcdM7MrMY5EZiZ1TgnAjOzGudEYGZW45wIzMxqnBOBmVmNy7LonJmZlUHWVVOdCMzMRrDhqJrqU0NmZiPQcFZNdSIwMxvRpnJk85O9qqYe+a4ngYayvYMTgZnZCNRTPvu22xazdfWkXlVTt945mdtu+6eylc/2GEGV89R7ZqPb7NkwZeneLJvR3Ktq6uzZ5XsPF52rYn0HkbbcUc+UCZ56z2y0KUfV1GJF53xEUMXyp5PMjQm6F7WzbEYzK1aMo6mp0tGZWbnkctDURGb/rz1GUIU8naSZlZMTQVXzdJJmNnSZjhFImgVcCeSAayLi833WHwX8OzAN+ExE/PNA+/QYwU61PvWemZWuImMEknLAVcApQCewTtLyiNiQt9nTwEeBd2UVx2iWy0HrLXWsWDGOjo5jaFjkq4bMbPCyHCyeDmyOiIcAJN0INAOvJIKI2AZsk3RGhnGMalkPIpnZ6JflGMEEYEvecmfaNmiS5ktqk9S2ffv2sgRnZmaJLBOBCrTt1oBERCyNiMaIaBw/fvwQwzIzs3xZJoJOoD5veSLweIbvZ2ZmuyHLRLAOOFzSZEl7AucAyzN8PzMz2w2ZDRZHRJekC4FVJJePXhsR90takK5fIukNQBswDnhZ0seBKRHxXFZxmZlZb5mWmIiI24Hb+7QtyXv+e5JTRmZmViGuNTQMXCHUzEYyJ4KMDcc0c2ZmQ+FaQxnLrxDaM83chs4XWbGi0pGZmSWcCDLiCqFmVi2cCDLnCqFmNrJ5hrKMuUKomY0EnqGsglwh1MxGOieCYeAKoWY2knmMwMysxjkRmJnVOCcCM7Ma50RgZlbjnAjMzGqcE4GZWY2rictHXf3TzKx/oz4RuPqnmVlxo/bUUE/Rt7q6JtZueqJX9c+1G5+grq7JRd/MzBjFiWCnqRzZ/GSv6p9HvutJoKGiUZmZjRSjNhG0tLQQEdx222K2rp7Uq/rn1jsnc9tt/0RE+IjAzGpepolA0ixJmyRtlnRxgfWS9LV0/X2SppU7htmzYcqEvVk2o5kfXzKNZTOamTJxL2bPLvc7mZlVp8wGiyXlgKuAU4BOYJ2k5RGxIW+z2cDh6WMG8K/pz7Jx9U8zs+KyvGpoOrA5Ih4CkHQj0AzkJ4Jm4NuRTIqwRtK+kg6MiN+VMxBX/zQz61+Wp4YmAFvyljvTtsFug6T5ktoktW3fvr3sgZqZ1bIsE4EKtPWdDq2UbYiIpRHRGBGN48ePL0twZmaWyDIRdAL1ecsTgcd3YxszM8tQlolgHXC4pMmS9gTOAZb32WY58L706qG3AM+We3zAzMyKy2ywOCK6JF0IrAJywLURcb+kBen6JcDtwOnAZuBF4ANZxWNmZoUpuWCnekjaDjwKHAA8WeFwRjL3T3Hun4G5j4qrtv45JCIKDrJWXSLoIaktIhorHcdI5f4pzv0zMPdRcaOpf0ZtiQkzMyuNE4GZWY2r5kSwtNIBjHDun+LcPwNzHxU3avqnascIzMysPKr5iMDMzMrAicDMrMZVXSIYaI6DaifpWknbJK3Pa9tf0o8k/Sb9uV/eukvSvtgk6bS89uMk/Tpd9zVJSttfJek7aft/S5qU95p56Xv8RtK8YfrIgyKpXtJPJD0g6X5JH0vb3UeApLGS1kq6N+2fy9N2908eSTlJ7ZJa0+Xa7p+IqJoHyR3KvwUOBfYE7gWmVDquMn/GE4BpwPq8ti8CF6fPLwa+kD6fkvbBq4DJad/k0nVrgeNJCvutAGan7R8BlqTPzwG+kz7fH3go/blf+ny/SvdHgf45EJiWPv8L4MG0H9xHSYwCXpM+HwP8N/AW988u/fQJ4HqgNV2u6f6peACD/Mc7HliVt3wJcEml48rgc06idyLYBByYPj8Q2FTo85OU8zg+3WZjXvu5wDfzt0mf15HcGan8bdJ13wTOrXRflNBXPyCZ/Mh9tGvf7AXcQzLZk/tnZ1wTgdXATHYmgprun2o7NVTS/AWj0OsjLcaX/nxd2t5ff0xIn/dt7/WaiOgCngVeW2RfI1Z6yD2V5Fuv+yiVnvboALYBP4oI909vXwU+Bbyc11bT/VNtiaCk+QtqSH/9Uayfduc1I46k1wA3AR+PiOeKbVqgbVT3UUR0R0QDyTff6ZKOKbJ5TfWPpCZgW0TcXepLCrSNuv6ptkRQq/MXPCHpQID057a0vb/+6Eyf923v9RpJdcA+wNNF9jXiSBpDkgSui4ib02b3UR8R8QfgLmAW7p8ebwPeKekR4EZgpqRl1Hr/VPrc1CDP7dWRDLBMZudg8dGVjiuDzzmJ3mMEX6L3QNYX0+dH03sg6yF2DmStIxkk7BnIOj1tv4DeA1nfTZ/vDzxMMoi1X/p8/0r3RYG+EfBt4Kt92t1HSYzjgX3T568G/gtocv8U7KsT2TlGUNP9U/EAduMf73SSK0V+C3ym0vFk8PluAH4H7CD5BvG3JOcXVwO/SX/un7f9Z9K+2ER61ULa3gisT9d9nZ13kY8FvkcyB8Ra4NC813wwbd8MfKDSfdFP/7yd5HD6PqAjfZzuPnolvjcB7Wn/rAcuS9vdP7v21YnsTAQ13T8uMWFmVuOqbYzAzMzKzInAzKzGORGYmdU4JwIzsxrnRGBmVuOcCKzsJJ0pKSQdleF7/HGIr78rrSbZkT7OLldslSbp1ZJ+Kik3hH28X9LX85YPlHRHke3vzK/YadXFicCycC7wc5KbaSpOiUK/6+dHREP6+H6f1+z2H9ER4IPAzRHRnd84xM80i6SYWn/+k6TqplUhJwIrq7QG0NtIboQ7J6/9xPRb+PclbZR0XV799tPTtp+ndd17asS3SPpk3j7W59d273k/Sasl3ZPWhm9O2ycpmbPgGyQVOPNv7e8v9kckXSbp58C7JZ0q6Vfpvr+XfraeOTEGFa+kuUrmCeiQ9M2eP8qS/ijpc0rmD1gj6fVp++sl3ZK23yvprZIWK51/Id3mc5I+WuCjnE9SlbWn338i6Xrg12nbrZLuVjJfwfy8/X1A0oOSfpr+G+abBaxIjwx+ln6O9ZL+Kl2/nOQLgFUhJwIrt3cBKyPiQeBpSdPy1k0FPk5S4/1Q4G2SxpKU450dEW8nKZEwGH8GzoyIacBJwJd7EgxwJPDtiJgaEY8WeO11eaeGXtuzvzSOO4FLgZPTfbcBn0jjvRqYA/wV8IaBApT0RuA9wNsiKQbXTfLHGmBvYE1EHAv8DPi7tP1rwE/T9mnA/cC/AfPSfe5Bkmiv6/Nee5LcyfpIXvN0krvwp6TLH4yI40jujP2opNem9XUuJ0kAp5D8G/XsMwccGREbgPNISsE3AMeS3NlNRDwDvCqvH62K1FU6ABt1ziUp8wtJUa9zSb6RA6yNiE4AJWWSJwF/BB6KiIfTbW4AXvmWWgIB/1fSCSRlhScAr0/XPRoRa4q89vyIaHtlR0n++E66+BaSP4a/SNv3BH4FHAU8HBG/SV+zrIR43wEcB6xL9/VqdhY1ewloTZ/fTfJHGJJa+e+DpJooSSnjZyU9JWlq+hnbI+KpPu91APCHPm1r8/oXkj/+Z6bP64HDSRLaXRGxPf1c3wGOSLeZQVLqG5L6OtcqKfx3a0R05O13G3AQ0DcmG+GcCKxs0m+DM4FjJAXJjHIh6VPpJv+Tt3k3ye9fodK8PbrofdQ6tsA255McRRwXETuUVJXs2e6FQX+Ina8RSS3/Xqc7JDXQf+ng/uIV8K2IuKTAa3bEzjovPX1SzDXA+0n+cF9bYP2f2LWfXukHSScCJ5NMnPKipLvytu/vc80GVgJExM/SpHsG8J+SvhQR3063G5u+v1UZnxqycjqb5FTMIRExKSLqSSosvr3IazYCh+ad+39P3rpHSE6LkJ5imlzg9fuQ1JffIekk4JChfYRXrCE5dfWX6fvvJemINN7Jkg5Lt8tPFP3Fuxo4W9Lr0nX7SxooztXAh9Ptc5LGpe23kJyvfzMFBm/TUzS59BRWIfsAz6RJ4CiSIx9IvvGfmJ4mGgO8O+8170jjIY17W0RcTXKqqufziiQ5PTLA57IRyInAyulckj9U+W4iOa9cUET8ieRqk5XpIO0TJKdBel67f3oa6cMkVWf7ug5olNRGcnSwcSgfIC+u7STfvG+QdB9JYjgqIv5Mciroh2m8+WMPBeNNz61fCtyR7utHJFMdFvMx4CRJvyY5ZXR0uq+XgJ+QlDbu7ue1d9B/8l0J1KVxLE4/F5HMytVCcvrrTtLTeZLGk4yb9Ez+cyLQIakdOAu4Mm0/jmSso2uAz2UjkKuPWsVJek1E/DH9VnkV8JuI+H+VjqsU6amWT0ZE0zC93x4kf6Tf3TNOUWCbqcAnIuK9ZXi/ucDEiPj8ANtdCSyPiNVDfU8bfh4jsJHg7yTNIxmQbSe5isj6kDSFZGD5lv6SAEBEtKeXjOaKHDWUJCKWlbjpeieB6uUjAjOzGucxAjOzGudEYGZW45wIzMxqnBOBmVmNcyIwM6tx/x8R2D8jGpZmcQAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Plot the data using errorbar(x,y,e).\n", "plt.errorbar(omega, Vratio, errVratio, fmt = 'ko', markersize = 5,\\\n", " linewidth = 1.8,\\\n", " markeredgecolor = 'b',\\\n", " markerfacecolor = (.49, 1, .63),\\\n", " capsize = 5)\n", "plt.xlabel('Angular Frequency (rad/s)')\n", "plt.ylabel('Voltage Ratio');" ] }, { "cell_type": "code", "execution_count": 59, "id": "swedish-latter", "metadata": {}, "outputs": [], "source": [ "# To do the actual fit, we will use the 'curve_fit()' function from the\n", "# 'SciPy' module. This way of fitting is very nice because we\n", "# use it for all types of fit models (linear, polynomial, line-in-parameter fits,\n", "# and nonlinear fit). It is capable of doing both unweighted\n", "# and weighted fits and it will return uncertainties in the fit parameters.\n", "\n", "# The first step is to define a function for the model that we will fit our\n", "# data to. In this case, the model is a Lorentzian.\n", "\n", "# w0 is the resonance frequency\n", "# A is the amplitude (height of the resonance peak)\n", "# width is the width of the resonance\n", "# x represents the independent variable along the horizontal axis\n", "# (angular frequency in this case)\n", "def LRCFunc(x, A, width, w0):\n", " y = A/np.sqrt(1+(x/width)**2*(1-(w0/x)**2)**2)\n", " return y" ] }, { "cell_type": "code", "execution_count": 60, "id": "thorough-johns", "metadata": {}, "outputs": [], "source": [ "# Here is the weighted fit of the nonlinear model to the data.\n", "# 'start' is used as initial guesses at the values of the best-fit parameters. \n", "start = (0.7, 0.5e5, 2e5)\n", "a_fit, cov = curve_fit(LRCFunc, omega, Vratio, sigma = errVratio,\\\n", " p0 = start, absolute_sigma = True)" ] }, { "cell_type": "code", "execution_count": 61, "id": "naval-allergy", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The best-fit parameters are: \n", " A ± ΔA = 0.7243152456809248 ± 0.003647934113877031 \n", " width ± Δwidth = 66816.88081227308 ± 623.0614968262379 \n", " w0 ± Δw0 = 213875.00919865974 ± 310.60920671532733\n" ] } ], "source": [ "# Extract the fit parameters and their uncertainties.\n", "print('The best-fit parameters are:\\\n", " \\n A \\u00B1 \\u0394A =', a_fit[0], '\\u00B1', np.sqrt(np.diag(cov))[0],\\\n", " '\\n width \\u00B1 \\u0394width =', a_fit[1], '\\u00B1', np.sqrt(np.diag(cov))[1],\\\n", " '\\n w0 \\u00B1 \\u0394w0 =', a_fit[2], '\\u00B1', np.sqrt(np.diag(cov))[2])" ] }, { "cell_type": "code", "execution_count": 62, "id": "christian-parade", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA5aklEQVR4nO3dd3hUZfbA8e9JASmioqg0KS4dqaEpShNCCQaCSHFdF11ZVkVQdxXRHwRQUdcGKrqogIoCshCqkEREOqEYegcFggXUFQwopJzfH3PBAUIIJMOdcj7PM09m7ty5c+Yl5Mx93/eeV1QVY4wxoSvM7QCMMca4yxKBMcaEOEsExhgT4iwRGGNMiLNEYIwxIS7C7QAu1DXXXKMVK1Z0OwxjjAkoa9eu/VFVS+X0XMAlgooVK7JmzRq3wzDGmIAiInvP9Zx1DRljTIizRGCMMSHOEoExxoQ4SwTGGBPiLBEYY0yIs0RgQlJ8fDwi4tzCEYlB5P+cn+GnnitZsiT/+te/+OWXX9wO2RifsURgQlJ8fDyqSmam0rnNEarJcP5FUaoyjKIk4JlZ3Yn//e8hXn55K1dddfWp5BAfH+9y9MYUrIC7jsCYgjRvHhxISWOjNiGSTJ5jCLVlNSWqH+LXbXu4QxNJKvIS32Q+RPFr+7F69UpKly7tdtjGFCg7IzAhLTUV2h2dQSSZAESSyU26i993/8BGbcJLDGb1b3WpGFGJQ4eiuPvuu8nKynI5amMKliUCE9Lq14ekYl3IcE6OM4hgeaGWdMo4PTnc/tsMTpzowsKFtxAREXvaOIJ1FZlA59NEICLtRWS7iOwSkUE5PP8vEVnn3DaJSJaIlPRlTMbAH4PFnTuHsy19O7VI4QmepxYp/HxiNwl6+6nk8DuFmBh+P1Wpc2ocIbr5j2RmKqpqicAEPPHVUpUiEg7sANoCacBqoJeqbjnH/p2BR1W1dW7HjYqKUqs1ZArS3Xf/hWnTjtKw4d9YvvwtIJGiTKMs5ehCIhPpzuUcZxP1iCSTDCKoF/EVLybcREyM29EbkzcislZVo3J6zpdnBI2BXaq6R1VPAJOB2Fz27wVM8mE8xpxl7969TJnyCQ89VJFlyzqgOgfVDI5k3sGrsxtQYsRTtL77T3SVOad1FXXKnMOXX/7ibvDGFBBfJoKywH6vx2nOtrOISFGgPTDtHM/3FZE1IrLm0KFDBR6oCS1ZWTBnDowYAY88kgSEM3DgwNP2CQ+HmBh45hno2ROSzxhHmEE0Gzd+dOmDN8YHfJkIJIdt5+qH6gwsU9Wfc3pSVceqapSqRpUqlWM5bWNy9ccFZOGUiJjF453X8uuQkWyb1YDCWVO54YaKZw385jaOcIA0kpIG2mCxCQq+HCNoBsSrarTz+CkAVR2Zw74JwFRV/eR8x7UxApMfc+bA0F7bWZle+1R/f1SRDTz3aY1c+/uzsjzXHKxbB/XqwY037qBmzWrEx8czdOjQSxW+MRfNrTGC1UAVEakkIoWAnsCsHIK7AmgBzPRhLMYAOV830OH3Waxbl/vrvLuKYmKgRo2qxMTEMGbMGI4fP+77wI3xIZ8lAlXNBB4GEoGtwKequllE+olIP69duwJJqnrUV7EYc1JO1w0kFetKvXoXfqyBAwdy8OBBJk+eXLBBGnOJ+axryFesa8hcjPj4eIYNGwaEUZSEU1NDZxDNAdI4Rlcgm6FDh+a5z19VqVmzJiVLlmTZsmW+DN+YfHOra8gYv3GyyJxqFoczOvNLqTeYWbUsr85uwJHMO1DNytPFYd5VS8PCIti27UaWL299VtVSG0A2gcTOCEzISUlJoWnTpkyYMIF77733oo6RlQVdo4+yY8F27iCRzyI7Ufm2G0lILEZ4eAEHbEwBsDMCY7xMnTqVyMhIYmNzu74xd6eqluIpTJea0ZC0lQeYN68AAzXmErFEYEKKqjJ16lTatWvHlVdeedHHyWn2UbtjCeedfWSMP7JEYELKqlWr2LdvH927d8/XcXKafTQrrMNFzT4yxm2WCExImTZtWr66hXK72nhv1jd07hxug8Um4NhgsQkpderUoVSpUixYsCDfx/K+2rh48V08+mg1xo59hwceeCD/gRpTwGyw2Bjg22+/ZePGjbRr165Ajud9tfGAATdSrVoVPvnkvFVSjPE7lghMyEhOTgYgOjq6wI8tIvTu3ZtFixaRlpZW4Mc3xpcsEZiQkZiYyHXXXUedOnV8cvzevXujqlZywgQcSwQmJGRnZ5OcnEzbtm0JC/PNr/2f/vQnGjdubN1DJuBYIjAhITU1lR9//NEn3UInZWVB3bpPk5raibff3k9Wls/eypgCZYnAhITExEQA2rZtW6DHPXPBm0XvluVfFOX1Bw9SImLWqfpDNp3U+DObPmpCQsuWLTl8+DCpqak+OX5OC940Kb6Z4ZOq2gL3xi/Y9FET0n799VeWLVvm026hHEtOHLWSEyYwWCIwQW/hwoVkZmYW2PUDOcmp5MScyE5WcsIEBEsEJuglJSVRtGhRbrnllgI/dm4lJ74+scdKTpiAYGMEJuhVqVKFqlWrMnfuXJ++j3fJiW+//Yy33+7Mxo3rqV27tk/f15i8sDECE7L27NnDrl27fDo+cJJ3yYmhQxsSFgaffvqpz9/XmPyyRGCCWlJSEuCbshK5ue6662jRogVTp04l0M66TejxaSIQkfYisl1EdonIoHPs01JE1onIZhFZ5Mt4TOhJSkrihhtuoGrVqpf8ve+66y62bdvGpk2bLvl7G3MhfJYIRCQceAvoANQEeolIzTP2uRIYA9yhqrWA/K0WYoyXjIwMFixYQHR0NCJyyd8/Li6OsLAw6x4yfs+XZwSNgV2qukdVTwCTgTNXA+kNTFfVfQCqetCH8ZgQk5KSwpEjR3w6bTQ31157La1ateLTTz+17iHj13yZCMoC+70epznbvFUFrhKRL0VkrYj8JacDiUhfEVkjImsOHTrko3BNsElKSiIsLIw2bdq4FsNdd93Fjh072LBhg2sxGHM+vkwEOZ2Ln/m1KAJoCHQCooH/E5GzOnNVdayqRqlqVKlSpQo+UhOUEhMTadKkCVdddZVrMXTt2pXw8HCmTp3qWgzGnI8vE0EaUN7rcTng2xz2ma+qR1X1R2AxUNeHMZkQ8fPPP7N69WrXuoVOKlWqlHUPGb/ny0SwGqgiIpVEpBDQE5h1xj4zgVtFJEJEigJNgK0+jMmEiM8//xxVveTTRnNy1113sXPnTtavX+92KMbkyGeJQFUzgYeBRDx/3D9V1c0i0k9E+jn7bAXmAxuAVcB7qmpz7Uy+JSYmcsUVV9CoUSO3QznVPWSzh4y/shITJmicLPHw1VfK6NF9aNHiGNOm+ccf3+joaHbv3s3OnTtdmcpqjJWYMEHrzIVhHu+8lvShL3D1T/2ZP/3PfrMwTPfu3dm9e7fP1kMwJj8sEZiAFh8fj6oye3YW1YtXYxNNeYnBbKIpVYtWZfbsLFTV9URg3UPGn1kiMEEhp4Vh2v82028Whrn66qu5/fbbmTJlis0eMn7HEoEJCjktDJNYrKtfLQzTu3dvvvnmG1auXOl2KMacxhKBCWi5LQyzPX2bXy0M06VLFy677DImTZrkdijGnMZmDZmgkZUF3buPY+bMvUyePIi4uCKEh7sd1em6d+/O4sWLOXDgABEREW6HY0KIzRoyISE8HHbvHkWLFkvo3t3/kgBAr169OHjwR0aO3MiIETBnjieBGeMmSwQmaHz33Xds2LDBL64mPtPJLqxu3bpTlAQmDsnm1yEjebzzWkpEzPKbaa4mNFkiMEEjOTkZwPX6QjnxnuZaMaLCadNcqxWv7jfTXE1oskRggkZiYiLXXnstdev6b93C1FTonDX/tGmu0UcT/GaaqwlNlghMUMjOziY5OZm2bdsSFua/v9b160Oyn09zNaHHf//HGHMB1q1bx6FDh/xyfAACa5qrCT2WCExQSExMBKBt27YuR5Kzk2MEqlkcybyDvv8uwr85Ssyj/+NI5h2o2hiBcY8lAhMUkpKSqFu3Ltdff73boZxXeDg8/nh1Klb8mC1b/u2X01xNaLFEYAJeeno6y5Yt89tuoZyICHfffTfJyckcOHDA7XBMiLNEYALewoULycjI8Mtpo7m59957yc7OZuLEiW6HYkKcJQIT8BITEylatCjNmzd3O5QLUqVKFZo3b86ECROsIqlxlSUCE/Dmz59P69atKVy4sNuhXLC//vWvbNu2jVWrVrkdiglhlghMQNu1axe7d++mffv2bodyUbp3706RIkWYMGGC26GYEGaJwAS0k9NGAzURlChRgm7dujFp0iR+//13t8MxIcqniUBE2ovIdhHZJSKDcni+pYgcFpF1zm2IL+MxwWf+/PnceOON3HjjjW6HctH++te/cvjwYWbOnOl2KCZE+SwRiEg48BbQAagJ9BKRmjnsukRV6zm34b6KxwSf48eP88UXXwTs2cBJrVq1onz58tY9ZFzjyzOCxsAuVd2jqieAyUCsD9/PhJilS5dy7NixgE8EYWFh3HvvvSQlJZGWluZ2OCYE+TIRlAX2ez1Oc7adqZmIrBeReSJSK6cDiUhfEVkjImsOHTrki1hNAJo/fz6FChWiZcuWboeSb/fddx+qyvvvv+92KCYE+TIRSA7bzpws/RVQQVXrAm8AM3I6kKqOVdUoVY0qVapUwUZpAlZiYiK33norxYsXdzuUfKtUqRLt2rXj3XffJTMz0+1wTIjxZSJIA8p7PS4HfOu9g6oeUdV05/5nQKSIXOPDmEyQOHDgABs3bgz4biFv/fr148CBA8ydO9ftUEyI8WUiWA1UEZFKIlII6AnM8t5BRK4XEXHuN3bi+cmHMZkgEejTRnMSExNDmTJleOedd9wOxYQYnyUCVc0EHgYSga3Ap6q6WUT6iUg/Z7c7gU0ish4YDfRUu9be5MH8+fMpW7YstWrlOKwUkCIiIvjb3/5GYmIiX3/9tdvhmBAiefm7KyJ3ALc5Dxep6myfRpWLqKgoXbNmjVtvb/xARkYGpUqV4s477+S9995zO5wCtX//fipWrMiTTz7J888/73Y4JoiIyFpVjcrpufOeEYjISGAAsMW5PeJsM8YVS5Ys4fDhw8TExLgdSoErX748MTExvP/++5w4ccLtcEyIyEvXUCegraqOU9VxQHtnmzGumD17NoULF/bb1cjyq1+/fhw8eJBp06a5HYoJEXkdI7jS6/4VPojDmDxRVWbPnk3r1q0pVqyY2+H4RHR0NFWqVGHUqFFuh2JCRF4SwUggVUQmiMgHwFrAOi+NK7Zv387u3bvp3Lmz26H4TFhYGAMGDCAlJYUVK1a4HY4JAedNBKo6CWgKTHduzVR1sq8DMyYns2d75ikE4/iAt3vvvZcrrriC119/3e1QTAg4ZyIQkerOzwZAaTwXiO0HyjjbjLlksrJgzhx4662rqFTpYcqUKX/+FwWw4sWL88ADDzBt2jT27dvndjgmyOV2RvCY8/OVHG4v+zguY07JyoKu0UcZ0nMbd+09RLG0vnSNPkpWltuR+Vb//v0BePPNN12OxAS7cyYCVe3r3O2gqq28b0DHSxOeCWXx8fGICBERMWxfsI2UozfxEoP5KqMB2xZsJyIiBhEhPj7e7VB94oYbbiAuLo53332X9PR0t8MxQSwvg8XL87jNmAIVHx+PqjJ8+Bzi5HMi8RRjiySTbpLMiBFzUNWgTQQAjz76KL/88otVJTU+ldsYwfUi0hAoIiL1RaSBc2sJFL1UARpTvz4kFe1CBhEAZBBBYrGu1KvnblyXQrNmzbj11lt5+eWX7QIz4zO5nRFE4xkLKAe8yh/jA48Bg30fmgl1J7uGOncOZ9vR7dQihSd4nlqksD19G507hwd119BJgwcPJi0tjYkTJ7odiglS5601JCLdVNVvLnG0WkOhqU+fvzFlyhGeeOIToqIi6NABwsPdjurSUFUaNmxIeno6W7duJTxUPrgpULnVGoo434tVdZqIdAJqAZd5bbf1hc0lkZmZyezZM+jaNZr4+PP+ygYdEWHw4MF0796dadOmcdddd7kdkgkyeSk69w7QA+iPZ9Wx7kAFH8dlzCmLFi3ip59+olu3bm6H4pq4uDiqV6/O888/j1VqNwUtL7OGblbVvwD/U9VhQDNOX3nMGJ+aPn06RYsWDapFaC5UWFgYgwYNYv369baCmSlweUkEvzk/j4lIGSADqOS7kIz5Q3Z2NgkJCXTo0IGiRUN7slrv3r2pXLkyQ4YMITs72+1wTBDJSyKYIyJXAv/Gs9j8N4DVGjKXxIoVK/juu++Ii4tzOxTXRUZGEh8fT2pqKtOmTTtVdmPECM/PYL/S2vhOnlYoO7WzSGE8A8aZqnrUZ1HlwmYNhZbHHnuMt956i0OHDlGiRAm3w3FdVlYWderUISsLqpRN4dtVB2h3dAZJxbpQtkk5EhKLhcxsKnNhLnqFMhEpKyJRzuLz4FmL4ElgZwHHaMxZsrOzmTZtGm3btrUk4BgxYgRbtmxh+/ZK7PhiOyvTazNSB7EyvXZIlN0wvpHblcUDgXXAG8BKEbkXzyL0RYCGlyI4E9qWL1/Ovn376NWrl9uh+I34+Hiys7MpXboTsSSFZNkNU/ByOyPoC1RT1WZAF+BdoJOqPqqq3+Xl4CLSXkS2i8guERmUy36NRCRLRO68kOBNcPvkk08oUqQIsbGxbofiV0SEvn0bMZN2IVl2wxS83BLB76r6M4Cq7gN2qOrKvB5YRMKBt4AOQE2gl4jUPMd+LwKJFxK4CW4ZGRlMnTqVO+64g+LFi7sdjt84WXZj2LAmpHEgZMtumIKV22Wa5URktNfja70fq+oj5zl2Y2CXqu4BEJHJQCyw5Yz9+gPTgEZ5jtoEvc8//5wff/zRuoXOEB8ff+qP/Jo1qTRqNITltzzIq4Ma0KFDA8LDbeqQuXC5JYJ/nfF47QUeuyyeFc1OSgOaeO8gImWBrkBrckkEItIXT1cVN9xwwwWGYQLRpEmTuPLKK0P6IrLziYqqT58+pZg4MZbq1bcQHv4nt0MyAeqciUBVP8jnsSWnw57x+HXgSVXNEslp91OxjAXGgmf6aD7jMn7u2LFjJCQk0KNHDwoXLux2OH7tueee49NPP+XJJ59k2jS/qQ1pAkxeLii7WGmcXoqiHPDtGftEAZNF5BvgTmCMiHTxYUwmAMyZM4f09HR69+7tdih+r3Tp0gwaNIjp06ezaNEit8MxAeqCLii7oAOLRAA7gDbAAWA10FtVN59j/wnAHFX9b27HtQvKgl9sbCyrV69m//79VnI5D44dO0aNGjW4/PLLSU1NJTIy0u2QjB+66AvK8kNVM4GH8cwG2gp8qqqbRaSfiPTz1fuawPb9998zd+5c7rnnHksCeVS0aFHeeOMNNm/ezGuvveZ2OCYA5WVhmqrA28B1qlpbROoAd6jqs5ciwDPZGUFwe/nll/nXv/7F1q1bqV69utvhBJQuXbqQlJTEli1bqFixotvhGD+T3zOCd4Gn8FQdRVU3AD0LLjxjPFSV8ePH07RpU0sCF2H06NGICI88cr6Z3cacLi+JoKiqrjpjW6YvgjGhbfXq1WzZsoU+ffq4HUpAuuGGGxg2bBizZ89mxowZbodjAkheEsGPInIjztRPpwxEnkpMGHMhxo8fT5EiRejRo4fboQSsAQMGUKdOHR588EF+/vlnt8MxASIvieAh4D9AdRE5AAwE/uHLoEzo+e2335g0aRLdunXjiiuucDucgBUZGcn48eM5dOgQAwYMcDscEyDOmwhUdY+q3g6UAqqranNV/cbnkZmQkpCQwOHDh61bqAA0aNCAp59+mokTJ1oXkcmTvMwaeiyHzYeBtaq6zhdB5cZmDQWnFi1akJaWxs6dOwkL8+V1jqHhxIkTNGnShG+//ZbNmzdzzTXXuB2ScVl+Zw1FAf3w1A4qi6fmT0vgXRF5oqCCNKFr06ZNLF68mH79+lkSKCCFChViwoQJ/O9//+PBBx/EVxeOmuCQl/91VwMNVPVxVX0cT2IoBdwG/NWHsZkQ8fbbb1O4cGHrFipgdevWZfjw4UydOpVx48a5HY7xY3lJBDcAJ7weZwAVVPU34LhPojIh49dff+Wjjz7irrvusu4LH3jiiSdo06YN/fv3Z+vWrW6HY/xUXhLBJ3iWqhwqIkOBZcAkESnG2WsLGHNBPv74Y3799VcefPBBt0MJSmFhYXz00UcUL16cnj178vvvv7sdkvFDeSo6JyJRwC14SksvVVXXRmttsDh4qCr16tUjLCyMr776itxKkZv8mTdvHh07duTBBx/krbfecjsc44J8F51z/vBPAqYDB0XEVocxFy0rC+bMgQce2MeGDeXp1+8hSwI+1qFDBx5//HHGjBnDxIkT3Q7H+JnzJgIRuUNEdgJfA4ucn/N8HZgJLifX2hUJp0TELB7vvJaS739CVYbxWL9rEbG1dn1t5MiRtGjRggceeIDU1FS3wzF+JC/XEazHs5Tk56paX0RaAb1Ute+lCPBM1jUU2ObMgaG9trMyvTaRZJJBBE2Kb2b4pKrExLgdXfA7ePAgDRs2JDw8nDVr1tgAfQjJb9dQhqr+BISJSJiqLgTqFWSAJnSkpkK7ozOIdOoWRpJJ9NEE1q1zN65Qce2115KQkMD3339Pjx49yMy0+pEmb4ngFxEpDiwGPhaRUVj1UXOBTnYNDRkSw3S9nQxnuewMIpimbfm//4uxrqFLJCoqinfeeYcvvviCRx55hMxMZc4cGDHCc8aWleV2hOZSy0vXUDHgNzxJ427gCmCiqrpS2tC6hgJbVhbUr7qZ3/ccJ06SSSrWlXJNypKQWAxbkOzSiY+PZ9iwYUAYRUmgHGWJJYmZtCONAxyjK5DN0KFDLTkHidy6hvKSCF5U1SfPt+1SsUQQ2I4fP06FCpUpU+Z+4uKGU68edOiAJQEXZGdn06LFvzm49HY20dTGbIJcfscI2uawrUP+QjKh6oMPPuCHH77lxRdv5ZlnICbGkoBbwsLCaN36MbqQbGM2Ie6ciUBE/iEiG4FqIrLB6/Y1sOHShWiCRUZGBiNHjqRJkybcfvvtbocT0k6O2Qwf3pUZtLUxmxB3zq4hEbkCuAoYCQzyeurXvI4PiEh7YBQQDrynqi+c8XwsMALIxjMAPVBVl+Z2TOsaClzjx4/nvvvuY86cOXTq1MntcAyeMZuu0UfZv2I/bY8lMJP2lG9WjsQlpexMLchc1BiBiJTM7aDnSwYiEg7swNO1lAasxnP9wRavfYoDR1VVRaQO8Kmq5rpquSWCwJSZmUmNGjUoUaIEa9assSuJ/UhWFsybB8nJh5gwYSAlSixj6dJFVKhQwe3QTAHKLRFE5PK6tTjrFOOpMeRNgcrned/GwC5V3eMEMRmIxatQnaqme+1fzOv9TJCZMmUKu3btYvr06ZYE/Ex4uGesJiamFPfd9wQtW7akTZs2LFq0iLJly7odnrkEzjlGoKqVVLWyc6t0xu18SQA8i9js93qc5mw7jYh0FZFtwFzgvpwOJCJ9RWSNiKw5dOhQHt7a+JOsrCyeffZZateuTWxsrNvhmFzUrVuXefPm8cMPP9CiRQv27t3rdkjmEshT0Tmn3tDLzi2vk8py+tp31jd+VU1wuoO64BkvOPtFqmNVNUpVo0qVKpXHtzf+4uOPP2bbtm0MGTLEViALAE2bNiU5OZkff/yR2267jd27d7sdkvGxvBSdewEYgKdLZwswQERG5uHYaUB5r8flgG/PtbOqLgZuFBErfhJEjh8/zpAhQ2jQoAHdunVzOxyTR02bNuWLL77g6NGj3HbbbWzbts3tkIwP5eXrWUegraqOU9VxQHsgL1M+VgNVRKSSiBQCegKzvHcQkT+J02EsIg2AQsBPF/IBjH/7z3/+w969e3nhhRfsbCDANGjQgC+//JLMzExatGiBTdIIXnn9n3ml1/0r8vICVc0EHgYSga14ZgRtFpF+ItLP2a0bsElE1gFvAT3UVtkOGr/++ivPPvssrVu3tusGAlTt2rVZvHgxRYsWpWXLlsybZxXog1Fus4ZOGgmkishCPP3+twFP5eXgqvoZ8NkZ297xuv8i8GKeozUB5dVXX+XQoUOMHDnSZgoFsGrVqrFixQo6duxI586d+c9//sP999/vdlimAOV2ZfGbInKzqk4CmuJZnWw60ExVJ1+qAE1g+v7773n55ZeJi4ujcePGbodj8un6669n0aJF3H777fztb39j6NChZGdnux2WKSC5dQ3tBF4RkW+AgcA+VZ2pqt9fisBMYBs8eDDHjx/nhRdeOP/OJiBcfvnlzJ49mz59+jB8+HC6d+9Oenr6+V9o/F5u1xGMUtVmQAvgZ2C8iGwVkSEiUvWSRWgCzurVqxk/fjyPPvooVapUcTscU4AiIyN5//33efXVV5kxYwbNmjVjz549bodl8um8ZahP21mkPjAOqKOqrlQisRIT/k1VueWWW9izZw87duygRIkSbodkfCQ5OZkePXogIkyZMsUmBPi5fJWhFpFIEeksIh/jWbR+B57ZPsac5ZNPPmHFihWMHDnSkkCQa9u2LatXr6Z06dJER0czbNgwsmx5s4CUW9G5tkAvPNcMrAImAzNU9eilC+9sdkbgv44cOUKNGjUoU6YMKSkpdt1AiEhPT+ehhx7iww8/pGXLlnz88ceUKVPG7bDMGS72jGAwsAKooaqdVfVjt5OA8W/PPPMM3333HW+++aYlgRBSvHhxPvjgAyZMmMCqVauoV68e8+fPBzyVTW09ZP93QWME/sDOCPzTqlWraNq0KQ899BBvvPGG2+EYl2zZsoVWrVpx8OBBbD1k/5KvNYv9jSUC/5ORkUFUVBQ//fQTW7ZssbGBEHfs2DEGDx7MqFG7qCbD2ahNbD1kP5DfNYuNydVrr73Ghg0beOONNywJGIoWLcrrr79Onz5vcIcm2nrIAcASgcmXHTt2EB8fT2xsLF27dnU7HOMHTq6HPH58f2bSztZDDgDWNWQuWmZmJs2bN2fHjh1s2rTJZoqY05xcDzkt5QDtjiYwg3bs1/106f1fXnvtZa699lq3Qwwp1jVkfOKll14iJSWFt956y5KAOUt4OCQkFmP4pKoUH/4kz02tycDBa5g6dTLVqlVj7Nixdt2Bn7AzAnNR1q9fT6NGjejatSuTJ0+26qImz7Zu3co//vEPFi1aRL169Xjttddo2bKl22EFPTsjMAXq+PHj3HPPPVx99dWMGTPGkoC5IDVq1GDhwoV88skn/PTTT7Rq1Yq4uDh27drldmghyxKBuWBPPPEEGzdu5L333uPqq692OxwTgESEXr16sX37dp599lmSkpKoWbMmjz/+OD/9ZIsUXmqWCMwFSUhIYPTo0QwYMIBOnfKyYqkx51akSBGefvppdu7cyT333MNrr71GpUqVGDp0KIcPH3Y7vJBhicDk2TfffMN9991HVFQUL730ktvhmCBSunRp3n//fTZu3Ei7du0YPnw4lSpV4oUXXuDoUats42uWCEyenDhxgp49e5Kdnc2UKVMoVKiQ2yGZIFSrVi3++9//snbtWpo1a8ZTTz1F5cqVefHFF+0MwYcsEZg8eeKJJ0hJSeH999+ncuXKbodjglyDBg2YO3cuy5Yto27dugwaNIgKFSowePBgp46RKUg+TQQi0l5EtovILhEZlMPzd4vIBue2XETq+jIec3E++OADRo0axcCBA7nzzjvdDseEkJtvvpmkpCRWr15N27ZteeGFF6hQoQIPP/wwX3/9tdvhBQ9V9ckNCAd2A5WBQsB6oOYZ+9wMXOXc7wCknO+4DRs2VHPprFy5UgsXLqxt2rTRjIwMt8MxIW7btm16//33a2RkpIaFhWmXLl30iy++0OzsbLdD83vAGj3H31VfnhE0Bnap6h5VPYFnYZvYM5LQclX9n/NwJVDOh/GYC/Tdd98RFxdHmTJlmDJlChEREW6HZEJctWrVeO+99/j6668ZNGgQS5YsoXXr1tSpU4d3332XY8eOuR1iQPJlIigL7Pd6nOZsO5f78SyFeRYR6Ssia0RkzaFDhwowRHMuv/32G3FxcRw+fJiZM2fa9QLGr5QtW5bnnnuO/fv3M27cOCIiIujbty/lypXjscceY8uWLW6HGFB8mQhyutw0x3oWItIKTyJ4MqfnVXWsqkapalSpUqUKMESTk6ysLO6++y+sXHkNnTuvYu/em2xlKeOXihQpQp8+ffjqq69YsmQJbdu25c0336RWrVrcfPPNjBs3jvT09NNeY6umnc2XiSANKO/1uBzw7Zk7iUgd4D0gVlXtkkI/8Oij/yQx4R5uKjSSilNmM7TXdrpGH7X/MMZviQjNmzdnypQpHDhwgFdeeYVffvmF+++/n9KlS9O3b19WrlxJZqbSNfooQ3tt59jQF+13+6RzDR7k9wZEAHuASvwxWFzrjH1uAHYBN+f1uDZY7DtDhw5VQKGTVmWNniBCFfQEEVqFtQqdFNChQ4e6Haox55Wdna3Lli3TPn36aGRkZMj/bpPLYLHPRv9UNVNEHgYS8cwgGqeqm0Wkn/P8O8AQ4GpgjFO4LFPPUR3P+F6NGjUAqFnzz3Te+jmR+sfKUt0kmWLD5/DMM25GaEzeiQg333wzN998M6+//jrTp09nxAiI3ZN02qppcZJM8VD/3T5XhvDXm50R+MasWbM0LCws5L81meCU17PdgQMHuh2qz5DLGYGtR2BITk4mJiaGevXqkZycTLFiJU6tLBV9NIHEYl0p16QsCYnFCA93O1pjLp73qmnRRxOYVySWzOsykGL3sGnTegCioqKIi4ujW7duVK1a1eWIC05u6xFYIghxS5YsITo6mipVqrBw4UJKliwJeP7DzJsH69ZBvXrQoQOWBExQONfv9s6dO0lISGD69OmkpKQAntpHsbGxdOzYkaZNmxIewP8JLBGYHC1fvpz27dtTpkwZFi9ebGvIGuNIS0tjxowZTJs2jSVLlpCVlUXJkiWJjo6mY8eOtG/fnmuuucbtMC+IJQJzli+//JKYmBhKly7NwoULKVfOLuo2Jie//PILycnJzJ07l3nz5nHw4EFEhCZNmtCxY0fatm1LVFSU3195b4nAnCYpKYnY2FgqV67M559/TunSpd0OyZiAkJ2dzVdffcVnn33G3LlzWb16NapKiRIlaNGiBW3atKFNmzbUqlXL75ZwtURgTpk9ezZ33nknNWrUIDk5GbtS25iLd+jQIRYuXMiCBQtYsGABu3fvBuC6666jdevWtGnThlatWlGpUiXXE4MlAgPAhx9+yP3330/9+vVJTEzkqquucjskY4LK3r17TyWFBQsW8MMPPwBQpkwZmjdvTvPmzbn11lu56aabLvnAsyWCEKeqjBw5kqeffpo2bdowffp0SpQo4XZYxgQ1VWXLli0sXryYpUuXsmTJEvbv99ThvPzyy7n55pu59dZbad68OY0aNaJo0aI+jccSQQjLysqif//+vP3229x9992MGzfOlpk0xiX79u07lRSWLl3Kpk2bAAgPD+emm26icePGNGnShMaNG1OjRo0CPWuwRBCijh07Ru/evZk5cyZPPPEEI0eOJCzMVic1xl/8/PPPLF++nJUrV7Jq1SpWrVp1am3mYsWKERUVRZMmTWjYsDFHj97G/v3X0KCBXNR1PZYIQtC+ffuIjY1l/fr1jBo1iv79+7sdkjHmPLKzs9m1axcpKSmnEkNq6gYiM6ZQjrLEksSc8A5QVnj6+a00bFiPKlWq5OnMIbdE4N8TX81FWbp0KXFxcRw/fpw5c+bQsWNHt0MyxuRBWFgYVatWpWrVquzevZtVq1YBnahEWTbRlEgyeS5rCLX2pfDnP08EehEZGUnDhg2pX78+9evXp169etSuXZsiRYrk+X3tjCDIvPfeezz44INUrFiRWbNmUb16dbdDMsbkw4gRcGzoi4zUQae2DZIXSH/wL0RFJbJu3TpSU1NZt24dR44cATwJpUqVKtx0002nbnFxcXZGEOyOHz/OY489xpgxY4iOjmbSpEk2PdSYABYfH8+wYcOATlRlGMOJIJJMMohgurZl51sPAHMZOnQoixYtIjs7m2+++YbU1FTWr1/Pxo0bSU1NZdq0aZzvC7+dEQSBPXv2cNddd7F27Vr++c9/MnLkSL+/3N0YkzdnVky90GrAR48eZfPmzTRp0sQGi4NRVhYMGbKCV175goiIjXz4YS/i4mLdDssYU8AKohqwzRoKQr/9doL61Taj+7PpQjJJRbtQvll5WzPAGJOj3BKBTSoPMPHx8YgIRYvGofuz2URTXuQpVh27iW0LthMREYOIEB8f73aoxpgAYR3JASQ7O5srr7ySwoULEx7enK6/2brCxpj8szOCAJGWlka7du149NFHOX78OMeOLSVBbyfDyeUZRDBN2/J//2dnBMaYC+PTRCAi7UVku4jsEpFBOTxfXURWiMhxEfmnL2MJVKrKBx98wE033cTKlSsZO3Ys2dnZZGbOoVqb6jQpvpmn5AWaFN9M9TbVyMycg6paIjDG5JnPuoZEJBx4C2gLpAGrRWSWqm7x2u1n4BGgi6/iCGRff/01f//730lOTqZ58+aMHz+eP/3pT4BnxkBCYjHmzavKunVPMryerStsjLk4vhwjaAzsUtU9ACIyGYgFTiUCVT0IHBSRTj6MI+BkZWUxevRonnnmGcLDwxkzZgx///vfzyoYFx4OMTGemzHGXCxfJoKywH6vx2lAk4s5kIj0BfoC3HDDDfmPzI9t3LiRBx54gJSUFGJiYhgzZgzly5d3OyxjTBDz5RhBTuuyXdRFC6o6VlWjVDUqWJdWPHz4MAMHDqR+/frs3r2bSZMmMWvWLEsCxhif82UiSAO8/4qVA7714fsFpOzsbD744AOqVq3K6NGj6du3L9u3b6dnz56ur3FqjAkNvuwaWg1UEZFKwAGgJ9Dbh+8XcNatW8dDDz3E8uXLadq0KfPmzaNBgwZuh2WMCTE+OyNQ1UzgYSAR2Ap8qqqbRaSfiPQDEJHrRSQNeAx4RkTSRCToF9P9/vvv6devHw0bNmTnzp2MGzeOZcuWWRIwxrjCp1cWq+pnwGdnbHvH6/73eLqMQsLRo0d59dVXefHFFzl+/Dj9+/cnPj6eK6+80u3QjDEhzEpM+FhWFsyZk8WECetZvHgUP/88kW7dujJy5EiqVKnidnjGGGOJwJcyM5UWjfbz84af6JydyLawgdRqOJIpU8rYhV/GGL9htYZ8ZPHixdStO4gf1x1iQ3ZjXmIwG7Ibk749nXnz3I7OGGP+YImggK1cuZIbb7yRFi1asGVLIWJJIpI/KoTenj6Nzp2fscJwxhi/YYmggKSmptK5c2eaNWvG999/f3IrM2l3WoXQGUQD69wK0xhjzmKJIJ/WrVvHnXfeSYMGDVi6dCnPP/88P/zwA6pqFUKNMQHBlqq8SCtWrOC5555j7ty5lChRggEDBvDYY4+dNRW0INYaNcaY/MptqUqbNXQBVJWFCxfy7LPPsnDhQq6++mpGjBjBww8/fM5rAaxCqDHG31kiyANVZe7cuTz33HOsXLmS0qVL88orr9C3b1+KFy/udnjGGJMvlghykZGRwdSpU/n3v//NunXrqFChAmPGjKFPnz5cdtllbodnjDEFwhJBDo4cOcK7777LqFGj2L9/P9WrV2fChAn07t2byMhIt8MzxpgCZYnAy/79+xk9ejRjx47lyJEjtGzZkjFjxtCxY8ezVgczxphgYYkAzxTQV155hcmTJ6OqdO/enccff5yoqBwH2I0xJqiEbCJQVRITE3n55ZdZsGABxYsX5+GHH2bAgAFUrFjR7fCMMeaSCfpEcHIef2oq1K8PzZv/yscff8gbb7zB9u3bKVOmDC+++CJ9+/a1ctDGmJAU1IkgKwu6Rh/lQEoabY/OYFBEB/Zlp/Fr1iM0atSQDz/8kB49elCoUCG3QzXGGNcEdSKYOzebfcv2svr3ukSSyYiMZ6gbsZb7X97C449Xczs8Y4zxC0E5FaZ58+aICLGxQ2n3+8zTqn/GZM7ln//8CBGhZcuW7gZqjDF+IKgSwY4dO3jkkUfYsGEDAFWrHiOx8B2nVf/8vHg3Zs9+FlXlyy+/dDFaY4zxDwGfCLKzs5k3bx4dOnSgWrVqvPPOO8TGxpKSksKWLa9QoXnF06p/lmtSlg4d3I7aGGP8h0+rj4pIe2AUEA68p6ovnPG8OM93BI4Bf1XVr3I75snqo0eOHGHChAm8+eab7Ny5k+uvv55//OMf9O3bl+uvv/7U/lb90xhjcq8+6rNEICLhwA6gLZAGrAZ6qeoWr306Av3xJIImwChVbZLbcWvXrq2tWrViwoQJpKen07RpUx555BG6detms3+MMeYc3CpD3RjYpap7nCAmA7HAFq99YoEP1ZONVorIlSJSWlW/O9dBN2/ezM6dO+nRowf9+/enUaNGPvwIxhgT/HyZCMoC+70ep+H51n++fcoCpyUCEekL9HUeHj9x4sSmjz76iI8++qhgIw4e1wA/uh2En7M2Oj9ro9wFWvtUONcTvkwEksO2M/uh8rIPqjoWGAsgImvOdXpjPKyNzs/a6PysjXIXTO3jy1lDaUB5r8flgG8vYh9jjDE+5MtEsBqoIiKVRKQQ0BOYdcY+s4C/iEdT4HBu4wPGGGMKns+6hlQ1U0QeBhLxTB8dp6qbRaSf8/w7wGd4ZgztwjN9tE8eDj3WRyEHE2uj87M2Oj9ro9wFTfv49DoCY4wx/i/gryw2xhiTP5YIjDEmxAVUIhCR9iKyXUR2icggt+MpaCIyTkQOisgmr20lRSRZRHY6P6/yeu4ppy22i0i01/aGIrLReW60U8oDESksIlOc7SkiUtHrNfc677FTRO69RB/5gohIeRFZKCJbRWSziAxwtlsbOUTkMhFZJSLrnTYa5my3NjqDiISLSKqIzHEeh24bqWpA3PAMOO8GKgOFgPVATbfjKuDPeBvQANjkte0lYJBzfxDwonO/ptMGhYFKTtuEO8+tAprhuU5jHtDB2f4g8I5zvycwxblfEtjj/LzKuX+V2+2RQ/uUBho49y/HU8KkprXRaW0kQHHnfiSQAjS1NsqxrR4DPgHmOI9Dto1c/8e4gH+0ZkCi1+OngKfcjssHn7MipyeC7UBp535pYHtOnx/P7Kxmzj7bvLb3Av7jvY9zPwLPVZHivY/z3H/w1IVyvT3O01Yz8dSysjbKuX2KAl/huaLf2uj0tikHLABa80ciCNk2CqSuoXOVowh216lzbYXz81pn+7nao6xz/8ztp71GVTOBw8DVuRzLbzmn2vXxfOO1NvLidHmsAw4CyapqbXS214EngGyvbSHbRoGUCPJUjiKEnKs9cmuni3mN3xGR4sA0YKCqHslt1xy2BX0bqWqWqtbD8623sYjUzmX3kGsjEYkBDqrq2ry+JIdtQdVGgZQIQrUcxQ8iUhrA+XnQ2X6u9khz7p+5/bTXiEgEcAXwcy7H8jsiEoknCXysqtOdzdZGOVDVX4AvgfZYG3m7BbhDRL4BJgOtRWQiodxGbvdNXUCfXgSegZVK/DFYXMvtuHzwOSty+hjBvzl9AOsl534tTh/A2sMfA1ir8QwQnhzA6uhsf4jTB7A+de6XBL7GM3h1lXO/pNttkUPbCPAh8PoZ262N/miLUsCVzv0iwBIgxtronO3Vkj/GCEK2jVz/h7jAf7SOeGaK7AaedjseH3y+SXhKcGfg+eZwP55+xQXATudnSa/9n3baYjvObAVnexSwyXnuTf64gvwyYCqekh6rgMper7nP2b4L6ON2W5yjfZrjOY3eAKxzbh2tjU5rozpAqtNGm4AhznZro5zbqyV/JIKQbSMrMWGMMSEukMYIjDHG+IAlAmOMCXGWCIwxJsRZIjDGmBBnicAYY0KcJQJT4ESkq4ioiFT34Xuk5/P1XzqVJNc5tzsLKja3iUgREVkkIuH5OMZfReRNr8elRSQpl/0/967WaQKLJQLjC72ApXgupHGdeOT0u363qtZzbv894zUX/UfUD9wHTFfVLO+N+fxM7fEUUjuXj/BU3DQByBKBKVBOHaBb8FwM19Nre0vnW/h/RWSbiHzsVbu9o7NtqVPT/WR9+HgR+afXMTZ513U/+X4iskBEvnLqwsc62yuKZ92CMXgqcHpf1n+u2L8RkSEishToLiLtRGSFc+ypzmc7uS7GBcUrIn8WzzoB60TkPyf/KItIuog8J571A1aKyHXO9utEJMHZvl5EbhaREeKsweDs85yIPJLDR7kbT2XWk+2+UEQ+ATY622aIyFrxrFfQ1+t4fURkh4gscv4NvbUH5jlnBoudz7FJRG51np+F5wuACUCWCExB6wLMV9UdwM8i0sDrufrAQDz13SsDt4jIZXhK8XZQ1eZ4SiRciN+BrqraAGgFvHIywQDVgA9Vtb6q7s3htR97dQ1dffJ4ThyfA88AtzvHXgM85sT7LtAZuBW4/nwBikgNoAdwi3qKwWXh+WMNUAxYqap1gcXAA8720cAiZ3sDYDPwPnCvc8wwPIn24zPeqxCeq1i/8drcGM+V+DWdx/epakM8V8U+IiJXO7V1huFJAG3x/BudPGY4UE1VtwC98ZSDrwfUxXN1N6r6P6CwVzuaABLhdgAm6PTCU+IXPAW9euH5Rg6wSlXTAMRTJrkikA7sUdWvnX0mAae+peaBAM+LyG14SgqXBa5znturqitzee3dqrrm1IE8+WOK87Apnj+Gy5zthYAVQHXga1Xd6bxmYh7ibQM0BFY7xyrCHwXNTgBznPtr8fwRBk+d/L+Ap5oonjLGh0XkJxGp73zGVFX96Yz3ugb45Yxtq7zaFzx//Ls698sDVfAktC9V9ZDzuaYAVZ19muAp9w2e2jrjxFP8b4aqrvM67kGgDHBmTMbPWSIwBcb5NtgaqC0iimdVORWRJ5xdjnvtnoXn9y+nsrwnZXL6WetlOexzN56ziIaqmiGeipIn9zt6wR/ij9cInlr+p3V3iEg9zl02+FzxCvCBqj6Vw2sy9I86LyfbJDfvAX/F84d7XA7P/8bZ7XSqHUSkJXA7nkVTjonIl177n+tzdQDmA6jqYifpdgI+EpF/q+qHzn6XOe9vAox1DZmCdCeerpgKqlpRVcvjqa7YPJfXbAMqe/X99/B67hs83SI4XUyVcnj9FXhqy2eISCugQv4+wikr8XRd/cl5/6IiUtWJt5KI3Ojs550ozhXvAuBOEbnWea6kiJwvzgXAP5z9w0WkhLM9AU9/fSNyGLx1umjCnS6snFwB/M9JAtXxnPmA5xt/S6ebKBLo7vWaNk48OHEfVNV38XRVnfy8gic5fXOez2X8kCUCU5B64flD5W0ann7lHKnqb3hmm8x3Bml/wNMNcvK1JZ1upH/gqTx7po+BKBFZg+fsYFt+PoBXXIfwfPOeJCIb8CSG6qr6O56uoLlOvN5jDznG6/StPwMkOcdKxrPMYW4GAK1EZCOeLqNazrFOAAvxlDXOOsdrkzh38p0PRDhxjHA+F+pZkSseT/fX5zjdeSJSCs+4yckFgFoC60QkFegGjHK2N8Qz1pF5ns9l/JBVHzWuE5HiqprufKt8C9ipqq+5HVdeOF0t/1TVmEv0fmF4/kh3PzlOkcM+9YHHVPWeAni/PwPlVPWF8+w3Cpilqgvy+57m0rMxAuMPHhCRe/EMyKbimUVkziAiNfEMLCecKwkAqGqqM2U0PJezhjxR1Yl53HWTJYHAZWcExhgT4myMwBhjQpwlAmOMCXGWCIwxJsRZIjDGmBBnicAYY0Lc/wO679w+MN/+TwAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Plot the data.\n", "plt.figure()\n", "plt.errorbar(omega, Vratio, errVratio, fmt = 'ko', markersize = 5,\\\n", " linewidth = 1.8,\\\n", " markeredgecolor = 'b',\\\n", " markerfacecolor = 'r',\\\n", " capsize = 5)\n", "plt.xlabel('Angular Frequency (rad/s)')\n", "plt.ylabel('Voltage Ratio')\n", "\n", "# Plot the best-fit line.\n", "xx = np.arange(1, 8e5, 100)\n", "plt.plot(xx, LRCFunc(xx, a_fit[0], a_fit[1], a_fit[2]), 'k-')\n", "plt.axis((0, 4.6e5, 0, 0.75));" ] }, { "cell_type": "code", "execution_count": 63, "id": "demonstrated-attack", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "48.14520250462778\n" ] } ], "source": [ "# Let's try to make our own estimates of the uncertainties in the best-fit parameters\n", "# using the methods discussed in the PHYS 232 lectures.\n", "\n", "# Let's give the best-fit parameters some recognizable names.\n", "AStar = a_fit[0]\n", "widthStar = a_fit[1]\n", "w0Star = a_fit[2]\n", "\n", "# The first step is to calculate the Chi-Squared value using the best-fit parameters\n", "# which should correspond to the minimum Chi-Squared value (ChiSqZero). The 'del w'\n", "# notation is used to clear any value previously assigned to w.\n", "ChiSq = 0\n", "for i in range(len(Vratio)):\n", " w = omega[i]\n", " ChiSq = ChiSq + ((Vratio[i] - LRCFunc(w, AStar, widthStar, w0Star))/errVratio[i])**2\n", "ChiSqZero = ChiSq\n", "print(ChiSqZero)\n", "del w" ] }, { "cell_type": "code", "execution_count": 64, "id": "welcome-dependence", "metadata": {}, "outputs": [], "source": [ "# The next step is to vary one of the parameters away from the value that minmizes\n", "# ChiSq to see how ChiSq depends on the paramter value. We expect this variation\n", "# to result in a quadratice dependence. As an example, we'll try finding the\n", "# uncertainty in the resonance frequency w0. Note, however, the same procedure\n", "# can be applied to each of the parameters (A, width, w0).\n", "\n", "# First, pick a suitable step size by which to vary w0. This may require some\n", "# trail and error.\n", "w0step = 5e2" ] }, { "cell_type": "code", "execution_count": 65, "id": "infectious-france", "metadata": {}, "outputs": [], "source": [ "# Now use nested loops to recalculate ChiSq at many values of w0 that bracket\n", "# the best-fit value a_fit[2]\n", "steps = 500\n", "ChiSqList = []\n", "w0List = []\n", "for i in range(steps):\n", " w0 = w0Star + (i - steps/2)*w0step \n", " w0List = w0List + [w0]\n", " ChiSq = 0\n", " for j in range(len(Vratio)):\n", " w = omega[j]\n", " ChiSq = ChiSq + ((Vratio[j] - LRCFunc(w, AStar, widthStar, w0))/errVratio[j])**2\n", " ChiSqList = ChiSqList + [ChiSq]\n", " del w" ] }, { "cell_type": "code", "execution_count": 66, "id": "extensive-russia", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAD4CAYAAADYU1DBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAdYUlEQVR4nO3df4wc533f8fendw4t25VMUieVvSNumYhISxlpbC1opimCIKxFxjVM/SEDVzQVkRAgqiqtkzRIRRiIUvsfK2mrRGilgpFUUaohimUciAgiOwSVwv2DIbX0L4piFJ7Nk3gWI15ARlFbQMkx3/6xz0pzy729/b07s58XsNjZZ2fmZrRHfe77PM/MKiIwMzPrt78z7AMwM7Px4MAxM7OBcOCYmdlAOHDMzGwgHDhmZjYQk8M+gF679dZbo1QqDfswzMxy5cyZM38REVP9/BmFC5xSqUSlUhn2YZiZ5Yqk1/v9M9ylZmZmA+HAMTOzgXDgmJnZQDhwzMxsIBw4ZmY2EA4cMzMbCAdOUiqVkLTi4et5zMx6Z83AkfSUpCuSXmnw3q9KCkm3ZtoOSJqX9JqkXZn2uySdTe89KkmpfZ2k51P7KUmlzDZ7JV1Ij71dn20Tr7/+OhGx4vH6632flm5mNjZaqXCeBnbXN0raDHwKeCPTtg2YA+5M2zwmaSK9/TiwH9iaHrV97gOuRcQdwCPAw2lfG4CHgE8C24GHJK1v7/TaU1/hmJlZ76wZOBHxDeBqg7ceAX4NyH6D2x7gcES8GxEXgXlgu6RNwM0RcTKq3/j2DHBPZptDafkosDNVP7uA4xFxNSKuAcdpEHy9VF/hmJlZ73Q0hiPps8APIuI7dW9NA5cyrxdT23Rarm9fsU1ELANvAxub7KvR8eyXVJFUWVpa6uSUgBvHcWptZmbWvbYDR9KHgC8Av97o7QZt0aS9021WNkYcjIhyRJSnpjq/91yjMRuP45iZ9UYnFc6PAFuA70haAGaAb0r6e1SrkM2ZdWeAN1P7TIN2sttImgRuodqFt9q++spdamZm/dF24ETE2Yi4LSJKEVGiGgyfiIg/B44Bc2nm2RaqkwNOR8Rl4B1JO9L4zH3AC2mXx4DaDLR7gZfSOM/XgbslrU+TBe5ObX0xOzsL0HDSgLvVzMy618q06OeAk8CPSlqUtG+1dSPiHHAEeBX4GvBARFxPb98PPEF1IsH3gBdT+5PARknzwK8AD6Z9XQW+BLycHl9MbX2xsLDw3nItfGpef/11h46ZWZdUtK6jcrkcnX4fzlpToYv238rMrEbSmYgo9/Nn+E4DGfWVjZmZ9Y4DJ6NZtxp4LMfMrBsOnDq1oMlOh27UZmZm7XHg1MlWOTXZoHGVY2bWGQfOGuqvyXGVY2bWGQdOA9nxG9/I08ysNxw4DSwsLDSdseZuNTOz9jlwVpEdy6m/3Y271czM2jc57APIA3epmZl1zxVOC3xNjplZ9xw4TTS7/sb3VzMza48Dp4lG1+RkeSzHzKx1DhwzMxsIB84a1rqhp7vVzMxa48BZw2rX5Pj+amZm7XHgtGCt+6uZmdnaHDht8hRpM7POOHBa5K8tMDPrjgOnRf7aAjOz7qwZOJKeknRF0iuZtt+S9KeSvivp9yV9NPPeAUnzkl6TtCvTfpeks+m9R5XuFyNpnaTnU/spSaXMNnslXUiPvb066X5wlWNm1lwrFc7TwO66tuPAxyLix4A/Aw4ASNoGzAF3pm0ekzSRtnkc2A9sTY/aPvcB1yLiDuAR4OG0rw3AQ8Ange3AQ5LWt3+KveMp0mZmnVszcCLiG8DVurY/iojl9PJPgJm0vAc4HBHvRsRFYB7YLmkTcHNEnIzqLZefAe7JbHMoLR8FdqbqZxdwPCKuRsQ1qiFXH3wD5SnSZmad68UYzi8AL6blaeBS5r3F1DadluvbV2yTQuxtYGOTfd1A0n5JFUmVpaWlrk5mLZ4ibWbWma4CR9IXgGXgK7WmBqtFk/ZOt1nZGHEwIsoRUZ6ammp+0H3mbjUzs8Y6Dpw0iP8Z4F9E7ZvJqlXI5sxqM8CbqX2mQfuKbSRNArdQ7cJbbV9D12wsx3eRNjNrrKPAkbQb+PfAZyPi/2XeOgbMpZlnW6hODjgdEZeBdyTtSOMz9wEvZLapzUC7F3gpBdjXgbslrU+TBe5ObUPnu0ibmbVvzW/8lPQc8NPArZIWqc4cOwCsA46n2c1/EhH/KiLOSToCvEq1q+2BiLiednU/1RlvN1Ed86mN+zwJPCtpnmplMwcQEVclfQl4Oa33xYhYMXnBzMzyQ+/3hhVDuVyOSqXS959TKpWaVjKzs7NrVkJmZqNC0pmIKPfzZ/hOAx1abYp0jcdyzMxWcuB0wWM5Zmatc+CYmdlAOHC65NvdmJm1xoHTpVbGcszMzIHTE2uN5bjKMTNz4PRFRJCdbu4qx8zMgdMz2W41SaQLYt/jKsfMxp0Dp0c8RdrMrDkHTg95xpqZ2eocOD3kL2gzM1udA6fH/AVtZmaNOXAGzN1qZjauHDh94Jt6mpndyIHTB56xZmZ2IwfOADSqeFzlmNm4ceD0STZkahVNozYzs3HhwOmTRlOkHTJmNs4cOH3km3qamb1vzcCR9JSkK5JeybRtkHRc0oX0vD7z3gFJ85Jek7Qr036XpLPpvUeVbjYmaZ2k51P7KUmlzDZ708+4IGlvz856gDxjzcysqpUK52lgd13bg8CJiNgKnEivkbQNmAPuTNs8JmkibfM4sB/Ymh61fe4DrkXEHcAjwMNpXxuAh4BPAtuBh7LBlheesWZmVrVm4ETEN4Crdc17gENp+RBwT6b9cES8GxEXgXlgu6RNwM0RcTKq9+1/pm6b2r6OAjtT9bMLOB4RVyPiGnCcG4MvdzxjzczGVadjOLdHxGWA9Hxbap8GLmXWW0xt02m5vn3FNhGxDLwNbGyyrxtI2i+pIqmytLTU4Sn1z1qz09y1ZmbjoNeTBtSgLZq0d7rNysaIgxFRjojy1NRUSwc6SGt9DTW4a83Miq/TwHkrdZORnq+k9kVgc2a9GeDN1D7ToH3FNpImgVuoduGttq9cWmssx8ys6DoNnGNAbdbYXuCFTPtcmnm2herkgNOp2+0dSTvS+Mx9ddvU9nUv8FIa5/k6cLek9WmywN2pLbf8fTlmNs4m11pB0nPATwO3SlqkOnPsy8ARSfuAN4DPAUTEOUlHgFeBZeCBiLiednU/1RlvNwEvpgfAk8CzkuapVjZzaV9XJX0JeDmt98WIqJ+8kCsLCwuUSqVVu8/crWZmRaZqMVEc5XI5KpXKsA+jqXQJUkOzs7PufjOzgZN0JiLK/fwZvtPAEPhiUDMbRw6cIfDFoGY2jhw4I8pVjpkVjQNnSHxdjpmNGwfOkLRyMairHDMrEgfOELUyluPQMbOicOAMmbvWzGxcOHCGrJVrblzlmFkROHBGgKscMxsHDpwR4AkEZjYOHDgjwhMIzKzoHDgjxF1rZlZkDpwR4gkEZlZkDpwR4yrHzIrKgTNiWplAIMmVjpnljgNnBLXSteZKx8zyxoEzotaqcsDjOWaWLw6cEdVqlePQMbO8cOCMsFaqHHetmVledBU4kn5Z0jlJr0h6TtIHJW2QdFzShfS8PrP+AUnzkl6TtCvTfpeks+m9RyUpta+T9HxqPyWp1M3x5k0rEwjAXWtmlg8dB46kaeDfAuWI+BgwAcwBDwInImIrcCK9RtK29P6dwG7gMUkTaXePA/uBremxO7XvA65FxB3AI8DDnR5vXrUSOq5yzCwPuu1SmwRukjQJfAh4E9gDHErvHwLuSct7gMMR8W5EXATmge2SNgE3R8TJiAjgmbptavs6CuysVT/jxPdaM7Mi6DhwIuIHwH8E3gAuA29HxB8Bt0fE5bTOZeC2tMk0cCmzi8XUNp2W69tXbBMRy8DbwMb6Y5G0X1JFUmVpaanTUxppvteameVdN11q66lWIFuAvw98WNLPNdukQVs0aW+2zcqGiIMRUY6I8tTUVPMDzzF3rZlZnnXTpfZPgYsRsRQRfwN8FfjHwFupm4z0fCWtvwhszmw/Q7ULbjEt17ev2CZ1290CXO3imHPN91ozszzrJnDeAHZI+lAaV9kJnAeOAXvTOnuBF9LyMWAuzTzbQnVywOnU7faOpB1pP/fVbVPb173AS2mcZ2y5yjGzvJrsdMOIOCXpKPBNYBn4FnAQ+AhwRNI+qqH0ubT+OUlHgFfT+g9ExPW0u/uBp4GbgBfTA+BJ4FlJ81Qrm7lOj7coFhYWKJVKTYOlVCq1VA2ZmQ2SilYwlMvlqFQqwz6Mvltrst7s7KxDx8xaJulMRJT7+TN8p4GccteameWNAyenPIHAzPLGgZNjrnLMLE8cODnmOxCYWZ44cHLOdyAws7xw4BSAu9bMLA8cOAXgCQRmlgcOnIJopcpx6JjZMDlwCsLfm2Nmo86BUyDuWjOzUebAKRhXOWY2qhw4BeNrc8xsVDlwCsjX5pjZKHLgFJS71sxs1DhwCsoTCMxs1DhwCqxRlZNtc5VjZoPkwCmwRhMIHDJmxVMqlZC04jGKPRgOnIJbq2ttFH8pzaw1taCp/0NyYmJiJP+4dOCMgWYTCDxjzSx/JicnGwZNzfXr1wd8RK3pKnAkfVTSUUl/Kum8pJ+QtEHScUkX0vP6zPoHJM1Lek3Srkz7XZLOpvcelaTUvk7S86n9lKRSN8c7rlqZJm1mo69W0YxqoKyl2wrnd4CvRcQ/AP4RcB54EDgREVuBE+k1krYBc8CdwG7gMUkTaT+PA/uBremxO7XvA65FxB3AI8DDXR6v0bjicZVjNpqy4zN5/+Ow48CRdDPwU8CTABHx1xHxl8Ae4FBa7RBwT1reAxyOiHcj4iIwD2yXtAm4OSJORkQAz9RtU9vXUWBnrfqx9qw1Oy3vv8hmRbPa+Eyr1roWbxi6qXB+GFgC/rukb0l6QtKHgdsj4jJAer4trT8NXMpsv5japtNyffuKbSJiGXgb2Fh/IJL2S6pIqiwtLXVxSsXlW96Y5UO3QQPVSQOtXIs3aN0EziTwCeDxiPg48H9J3WeraFSZRJP2ZtusbIg4GBHliChPTU01P+ox5lvemI2utSYCtGp2dpbl5eUeHVVvdRM4i8BiRJxKr49SDaC3UjcZ6flKZv3Nme1ngDdT+0yD9hXbSJoEbgGudnHMY8+3vDEbHdnxmW4mAkxMTBARRMRIVjY1HQdORPw5cEnSj6amncCrwDFgb2rbC7yQlo8Bc2nm2RaqkwNOp263dyTtSOMz99VtU9vXvcBLaZzHOuRb3piNhsnJya7/wKsFzahWNPUmu9z+3wBfkfRDwPeBn6caYkck7QPeAD4HEBHnJB2hGkrLwAMRUYv0+4GngZuAF9MDqhMSnpU0T7WymevyeI1qldPsF91Vjll/lEqlnvz7mp2dHelKZjUqWsFQLpejUqkM+zBGXqNf/GwQFe33wmyYehE0ExMTfa1kJJ2JiHLffgC+08DYavTXUfYfhLvVzLo3DhMB2uHAGWO+5Y1Zf/TijgB5mQjQjm7HcCzHFhYWaHYdrcdyzNozOTnZ9W1n8jo+0wpXOGPOF4OadadXU5tnZ2cLVc004sAZc/W/3I2+P8ehY3ajXt0RoGjdZs04cMz3WTNrgycCdM6BY74Y1GwN43ZHgH5x4BjgW96YrWYc7wjQLw4cA3w3abOsXk8EGPegqXHg2Ht8N2kbd54I0F8OHFvBXWs2jnoRNDCeEwHa4Qs/bYW1LgaF6j9O/+VmRdHtxZr9vsdZkbjCsRu4yrFx0O3tZzw+0z4Hjt3AEwis6LqZeTYOdwToFweONeQJBFZE3VQ1DpruOXBsVe5asyLptKpx0PSOA8dW5TsQWBF0WtU4aHrPgWNNucqxPOvkmzYdNP3jwLGmWplAIMmVjo2cTrrQivxdNKOg68CRNCHpW5L+IL3eIOm4pAvpeX1m3QOS5iW9JmlXpv0uSWfTe48qXQgiaZ2k51P7KUmlbo/X2tfKP0BXOjZK2r22xlXNYPSiwvk8cD7z+kHgRERsBU6k10jaBswBdwK7gcckTaRtHgf2A1vTY3dq3wdci4g7gEeAh3twvNaBtaoc8HiODV8n4zWuaganq8CRNAP8M+CJTPMe4FBaPgTck2k/HBHvRsRFYB7YLmkTcHNEnIyIAJ6p26a2r6PAzlr1Y4PVapXj0LFhaXe8xlXN4HVb4fw28GvA32babo+IywDp+bbUPg1cyqy3mNqm03J9+4ptImIZeBvYWH8QkvZLqkiqLC0tdXlKtppWqhx3rdkwdBI2DprB6zhwJH0GuBIRZ1rdpEFbNGlvts3KhoiDEVGOiPLU1FSLh2PtamUCAbhrzQarnckBtTs5O2yGo5sK5yeBz0paAA4DPyPpfwBvpW4y0vOVtP4isDmz/QzwZmqfadC+YhtJk8AtwNUujtm61ErouGvNBqWdyQG+yebwdRw4EXEgImYiokR1MsBLEfFzwDFgb1ptL/BCWj4GzKWZZ1uoTg44nbrd3pG0I43P3Fe3TW1f96afcUOFY4Pl0LFha3dygL82YDT04+sJvgwckbQPeAP4HEBEnJN0BHgVWAYeiIjab8v9wNPATcCL6QHwJPCspHmqlc1cH47XOtDK1xh4PMf6xeM1+aSiFQzlcjkqlcqwD2MstDJQ63/s1kur/c5NTEw0rHb8+9c6SWciotzPn+E7DVjHPFXaBqnZHzj1YePJAaPJgWNdaXWqtEPHutHOtGdPDhhdDhzrysLCAq10y3o8xzrVzrRnTw4YbQ4c6wlfn2P9UCqV2pqJ5i600ebAsZ7weI71WjvdaA6bfHDgWM94PMd6pd1uNIdNPjhwrGfaGc9x6Nhq2rl7gMMmXxw41nONKp2JiYkVrx06Vq+duwd42nM+OXCs5xr9T6DR/0Q8c82y2vl98Ey0fHLgWF+0Mp4D1e4Ts3Z+D1r93bLR48Cxvmh1POf69evuWhtj7kYbLw4c66tWZ6650hk/vnvA+HHgWF+50rFGfPeA8eTAsYHwNTpW42nP48uBYwOxsLBww9ToRhw6xdbOrWomJiYcNgXjwLGBWV5edqUzxtyNZg4cG6h2Kh0rDnejGThwbAha/ctVkiudnPO0Z8ty4NhQtHrxnrvX8svTnq1ex4EjabOkP5Z0XtI5SZ9P7RskHZd0IT2vz2xzQNK8pNck7cq03yXpbHrvUUlK7eskPZ/aT0kqdXGuNkJanS4N7l7Lo3a/WsBhMx66qXCWgX8XEf8Q2AE8IGkb8CBwIiK2AifSa9J7c8CdwG7gMUm1zvzHgf3A1vTYndr3Adci4g7gEeDhLo7XRlCrlY671/KjnckBnok2XjoOnIi4HBHfTMvvAOeBaWAPcCitdgi4Jy3vAQ5HxLsRcRGYB7ZL2gTcHBEno/on7zN129T2dRTYWat+rBgWFhbcvVYQ7YzXgLvRxlFPxnBSV9fHgVPA7RFxGaqhBNyWVpsGLmU2W0xt02m5vn3FNhGxDLwNbGzw8/dLqkiqLC0t9eKUbIAcOvnXThcauBttXHUdOJI+Avwe8EsR8VfNVm3QFk3am22zsiHiYESUI6I8NTW11iHbCHLo5FcnYeNutPHUVeBI+gDVsPlKRHw1Nb+VuslIz1dS+yKwObP5DPBmap9p0L5iG0mTwC3A1W6O2UaXQyd/2hmvAYfNuOtmlpqAJ4HzEfGfM28dA/am5b3AC5n2uTTzbAvVyQGnU7fbO5J2pH3eV7dNbV/3Ai9Fq1ObLJfaDR3fZXp42rmY09fYGHRX4fwk8C+Bn5H07fT4NPBl4FOSLgCfSq+JiHPAEeBV4GvAAxFR+229H3iC6kSC7wEvpvYngY2S5oFfIc14s2JrJ3SuX7/u0BmwdicHeLzGalS0gqFcLkelUhn2YViPtDop0V01g+HxmuKSdCYiyv38Gb7TgI00d6+NDoeNdcv/Qm2kLSwstPw/ulr3mrtveq/d8Rp/BtaIKxwbebXb4LRyl+nr16/7rgQ95Is5rZccOJYby8vLLYUOeNp0L3Qy5dlhY804cCxXWv0SN/C4TqfarWrA4zXWGgeO5U6rX+IG7mJrV7tVDThsrHUOHMuldrrXwNXOWjqpanwxp7XLgWO51U73Gvgi0dV0UtV4coB1woFjudbOXQng/S42B081aDqtahw21gkHjuVeO9Oma8a52umk+ww8C82658Cxwmh3XGfcqp1a0HTyld2eGGC94MCxQlleXna1U6fboPHEAOsVB44Vkqud7oIGXNVY7zlwrLDancUGxQieXgSNqxrrBweOFVonEwognxeMdhs0vq7G+s2BY2Ohk2oHqheMjnrFU5ve3G3QeAaa9ZsDx8ZGp9UOvF/xjErVUwuZTqY31zhobNAcODZ2Oq12ampVz6DDpxchU+NramwYHDg2lmrVTjfBA/0Nn9qYTK9CBjwhwIZrdDumMyTtBn4HmACeiIgvD/mQrCBq/+Nt5xstV1MLn6xW7jnWi5/djO97ZqNi5CscSRPAfwV+FtgG/HNJ24Z7VFY0nVww2ors2M9qj36GjbvObJSMfOAA24H5iPh+RPw1cBjYM+RjsoKqBU+3XW3DVJsM4K4zGzV5CJxp4FLm9WJqe4+k/ZIqkipLS0sDPTgrpl6N8QxKbWzGs85slOUhcNSgLVa8iDgYEeWIKE9NTQ3osGwc1IKn9uh1l1s3siHjSsbyIA+TBhaBzZnXM8CbQzoWG3PZ6qHfg/31PPhveZeHCudlYKukLZJ+CJgDjg35mMzeG+/JPnrRBZcdg8k+HDaWdyNf4UTEsqRfBL5OdVr0UxFxbsiHZdaQu7bMVjfygQMQEX8I/OGwj8PMzDqXhy41MzMrAAeOmZkNhAPHzMwGwoFjZmYDoYhYe60ckbQEdPZNVJ27FfiLAf/MYfM5jwef83i4FfhwRPT1yvnCBc4wSKpERHnYxzFIPufx4HMeD4M6Z3epmZnZQDhwzMxsIBw4vXFw2AcwBD7n8eBzHg8DOWeP4ZiZ2UC4wjEzs4Fw4JiZ2UCMdeBIekrSFUmvZNo2SDou6UJ6Xp9574CkeUmvSdqVab9L0tn03qOSlNrXSXo+tZ+SVMpsszf9jAuS9g7olFc759+Q9ANJ306PT2feK8I5b5b0x5LOSzon6fOpvbCfdZNzLuxnLemDkk5L+k465/+Q2ov8Oa92zqP5OTf63o1xeQA/BXwCeCXT9pvAg2n5QeDhtLwN+A6wDtgCfA+YSO+dBn6C6reTvgj8bGr/18B/S8tzwPNpeQPw/fS8Pi2vH+I5/wbwqw3WLco5bwI+kZb/LvBn6dwK+1k3OefCftbp+D6Slj8AnAJ2FPxzXu2cR/JzHusKJyK+AVyta94DHErLh4B7Mu2HI+LdiLgIzAPbJW0Cbo6Ik1H9FJ6p26a2r6PAzvRXwy7geERcjYhrwHFgd6/Pr5FVznk1RTnnyxHxzbT8DnAemKbAn3WTc15NEc45IuL/pJcfSI+g2J/zaue8mqGe81gHzipuj4jLUP1HC9yW2qeBS5n1FlPbdFqub1+xTUQsA28DG5vsa5h+UdJ3Ve1yq3U5FO6cU3fAx6n+JTgWn3XdOUOBP2tJE5K+DVyh+j/Dwn/Oq5wzjODn7MBpnRq0RZP2TrcZhseBHwF+HLgM/KfUXqhzlvQR4PeAX4qIv2q2aoO2XJ53g3Mu9GcdEdcj4seBGap/uX+syepFPueR/JwdODd6K5WXpOcrqX0R2JxZbwZ4M7XPNGhfsY2kSeAWqt1Zq+1rKCLirfRL+7fA7wLb01uFOWdJH6D6P96vRMRXU3OhP+tG5zwOnzVARPwl8L+odvEU+nOuyZ7zyH7OvRzAyuMDKLFyAP23WDnA+Jtp+U5WDrZ9n/cH216mOlBXG2z7dGp/gJWDbUfi/cG2i1QH2tan5Q1DPOdNmeVfptrHW5hzTsf4DPDbde2F/aybnHNhP2tgCvhoWr4J+N/AZwr+Oa92ziP5Off9H/soP4DnqJabf0M1rfdR7Zs8AVxIzxsy63+B6qyO10gzOFJ7GXglvfdfeP8ODh8E/ifVgbnTwA9ntvmF1D4P/PyQz/lZ4CzwXeBY3S9rEc75n1At9b8LfDs9Pl3kz7rJORf2swZ+DPhWOrdXgF9P7UX+nFc755H8nH1rGzMzGwiP4ZiZ2UA4cMzMbCAcOGZmNhAOHDMzGwgHjpmZDYQDx8zMBsKBY2ZmA/H/AbBYiGsg7mS8AAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Plot the ChiSq value as a function of w0.\n", "plt.plot(w0List, ChiSqList, 'ks', linewidth = 2.5, markersize = 5, fillstyle = 'none');" ] }, { "cell_type": "code", "execution_count": 67, "id": "signed-conspiracy", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAASUUlEQVR4nO3df4xdaV3H8feX7m4XC0orbVPa2takCi2/GQsG/1BWt+VHbA3ZpCDayCb9w5pAokILJsIfm/DDKDG6wUbRIaKlApttiAKlgr8CW6bLQrddaoftdHds0xZWhZWka8vXP+4z9sx0pnPvzNyZe595v5Kbe+5zz7nnOzP3fu6Z55zznMhMJEl1edZCFyBJmnuGuyRVyHCXpAoZ7pJUIcNdkipkuEtShdoK94gYiYiTEfFIRAyVthURcTQizpb75Y35D0TEcESciYjt3SpekjS5TrbcfyEzX56ZA+XxfuBYZm4GjpXHRMQWYDewFdgB3B8RS+awZknSNGbTLbMTGCzTg8CuRvuhzLyameeAYWDbLNYjSerQbW3Ol8AXIiKBP8vMg8DqzLwIkJkXI2JVmXct8NXGsqOlbZyI2AvsBVi2bNmrXvjCF87wR5jcyZMneeaZZ8a13XHHHbzkJS+Z0/VI0kI5ceLEdzJz5WTPtRvur83MCyXAj0bEt24xb0zSdtMYB+UL4iDAwMBADg0NtVmKJAkgIs5P9Vxb3TKZeaHcXwYeoNXNciki1pQVrAEul9lHgfWNxdcBFzovW5I0U9OGe0Qsi4jnjk0DdwOPAkeAPWW2PcCDZfoIsDsilkbEJmAzcHyuC5ckTa2dbpnVwAMRMTb/32Tm5yLia8DhiLgXeAK4ByAzT0XEYeA0cA3Yl5nXu1K9JGlS04Z7Zj4OvGyS9u8Cd02xzH3AfbOuTpI0I56hKkldsnHjRiJi3G3jxo3zsu52j5aRJHXo/PnzTLwgUuni7jq33CWpQoa7JFXIcJekCtnnLkldsmHDhpv62Dds2DAv6zbcJalLRkZGFmzddstIUoUMd0mqkOEuSRUy3CWpQoa7JFXIcJekChnuklQhw12SKmS4S1KFDHdJqpDhLkkVMtwlqUKGOwt7KSxJ6gZHhWRhL4UlSd3glrskVchwl6QKGe6SVCH73FnYS2FJUjcY7izspbAkqRvslpGkChnuklQhw12SOtAvJz3a5y5JHeiXkx7dcpekChnuklQhw12SKmSfuyR1oF9OejTcJakD/XLSY9vdMhGxJCK+HhGfLY9XRMTRiDhb7pc35j0QEcMRcSYitnejcEnS1Drpc38H8Fjj8X7gWGZuBo6Vx0TEFmA3sBXYAdwfEUvmplxJUjvaCveIWAe8EfjzRvNOYLBMDwK7Gu2HMvNqZp4DhoFtc1KtJKkt7W65fwR4F/DDRtvqzLwIUO5Xlfa1wJON+UZL2zgRsTcihiJi6MqVK53WLUm6hWnDPSLeBFzOzBNtvuZkp2rlTQ2ZBzNzIDMHVq5c2eZLS5La0c7RMq8Ffjki3gDcCfxoRPw1cCki1mTmxYhYA1wu848C6xvLrwMuzGXRkqRbm3bLPTMPZOa6zNxIa0fpP2bm24AjwJ4y2x7gwTJ9BNgdEUsjYhOwGTg+55VLkqY0m+PcPwAcjoh7gSeAewAy81REHAZOA9eAfZl5fdaVSpLaFhNHN1sIAwMDOTQ0tNBlSFJfiYgTmTkw2XOOLSNJFTLcJalChrskVchw71C/XGJL0uLmqJAd6pdLbEla3Nxyl6QKGe6SVCHDXZIqZJ97h/rlEluSFje33Ds0MjJCZo679ctltyTdrNYj4Nxyl7So1XoEnFvuklQhw12SKmS4S1KF7HOXtKjVegSc4S5pUav1aDe7ZSSpQoa7JFXIcJekChnuklQhw12SKmS4S1KFDHdJqpDhLkkVMtwlqUKGuyRVyHCXpAoZ7l1U6xVeJPU+Bw7rolqv8CKp97nlLkkVMtwlqUKGuyRVyHDvorErvDRvNVzhRep1HszgDtWuqvUKL1Kv82CGNrbcI+LOiDgeEd+IiFMR8f7SviIijkbE2XK/vLHMgYgYjogzEbG9mz+AJOlm7XTLXAVel5kvA14O7IiI1wD7gWOZuRk4Vh4TEVuA3cBWYAdwf0Qs6ULtkqQpTBvu2fJ0eXh7uSWwExgs7YPArjK9EziUmVcz8xwwDGyby6IlSbfW1g7ViFgSEY8Al4GjmfkQsDozLwKU+1Vl9rXAk43FR0vbxNfcGxFDETF05cqVWfwIkjSeBzO0uUM1M68DL4+I5wEPRMSLbzH7ZHst8qaGzIPAQYCBgYGbnpekmfJghg4PhczM/wK+TKsv/VJErAEo95fLbKPA+sZi64ALsy1UktS+do6WWVm22ImIZwO/CHwLOALsKbPtAR4s00eA3RGxNCI2AZuB43NctyTpFtrpllkDDJYjXp4FHM7Mz0bEV4DDEXEv8ARwD0BmnoqIw8Bp4Bqwr3TrSJLmSUw80H8hDAwM5NDQ0EKXIUl9JSJOZObAZM85/IAkVchwl6QKGe6SVCHDXZIqZLhLUoUMd0mqkOEuSRUy3HuEV46RNJe8ElOP8MoxkuaSW+6SVCHDXZIqZLhL6gvul+qMfe49YuzKMRPbJLW4X6ozhnuP8MoxkuaS3TKSVCHDXZIqZLeMpL7gfqnOGO6S+oL7pTpjt4wkVchwl6QKGe6SVCHDXZIqZLhLUoUMd0mqkOEuSRUy3CWpQoa7JFXIcO9DjmstaToOP9CHHNda0nTccpekChnukhaMXYzdY7eMpAVjF2P3GO59yHGtJU3HcO9DjmstaTr2uUtShaYN94hYHxFfiojHIuJURLyjtK+IiKMRcbbcL28scyAihiPiTERs7+YPIKl/jXUxNm92Mc6NdrbcrwG/nZkvAl4D7IuILcB+4FhmbgaOlceU53YDW4EdwP0RsaQbxUvqbyMjI2TmuJvdjnNj2nDPzIuZ+XCZ/j7wGLAW2AkMltkGgV1leidwKDOvZuY5YBjYNsd1S5JuoaM+94jYCLwCeAhYnZkXofUFAKwqs60FnmwsNlraJr7W3ogYioihK1euzKB0SdJU2g73iHgO8GngnZn5vVvNOklb3tSQeTAzBzJzYOXKle2WIUlqQ1vhHhG30wr2T2TmZ0rzpYhYU55fA1wu7aPA+sbi64ALc1OuJKkd7RwtE8BfAI9l5h82njoC7CnTe4AHG+27I2JpRGwCNgPH565kSdJ02jmJ6bXArwEnI+KR0vYe4APA4Yi4F3gCuAcgM09FxGHgNK0jbfZl5vW5LlySNLVpwz0z/5XJ+9EB7ppimfuA+2ZRlyRpFjxDVZIqZLgvAg6rKi0+hvsiMDasavN2/vz5hS5LlXJjojc4KqSkOeUY7b3BLXdJqpDhLkkVsltmEfDKTdLiY7gvAg6hqvnkxkRvMNwlzSk3JnqDfe6SVCHDXZIqZLhLUoUMd0mqkOEuSRUy3CWpQoa7JFXIcNc4jugn1cFw1zgOD6zJ+KXffzxDVdK0HMa3/7jlLkkVMtwlqUJ2y2gcR/ST6mC4axxH9NNk/NLvP4a7pGn5pd9/7HOXpAoZ7pJUIcNdkipkuEtShQx3SaqQ4a4Zc7yR/ubfr24eCqkZc7yR/ubfr25uuUtShQx3SaqQ4S5JFZo23CPiYxFxOSIebbStiIijEXG23C9vPHcgIoYj4kxEbO9W4Vp4Y+ONNG+ON9I//PvVrZ0t978Cdkxo2w8cy8zNwLHymIjYAuwGtpZl7o+IJXNWrXrKyMjITVdtcgyS/uHfr27Thntm/jPw1ITmncBgmR4EdjXaD2Xm1cw8BwwD2+amVElSu2ba5746My8ClPtVpX0t8GRjvtHSdpOI2BsRQxExdOXKlRmWIUmazFzvUJ3sINmcpI3MPJiZA5k5sHLlyjkuQ5IWt5mG+6WIWANQ7i+X9lFgfWO+dcCFmZenWng2pDS/ZhruR4A9ZXoP8GCjfXdELI2ITcBm4PjsSlQNxs6GbN7Onz+/0GVVxS9QNU07/EBE/C3w88DzI2IU+H3gA8DhiLgXeAK4ByAzT0XEYeA0cA3Yl5nXu1S7pAaHE1DTtOGemW+Z4qm7ppj/PuC+2RQlSZodz1CVpAo5KqTmxdjZkBPbJHWH4a554ZmP3ecXqJoMd6kSfoGqyT539RwP6ZNmz3BXz/GY+Bv8otNM2S0j9TCPXddMueUuSRUy3CWpQnbLqOd4SJ80e265q+d0coWgftzh2EnNXgpPM+WWu/paP+5w7KRmj13XTLnlLkkVMtwlqUKGuxaVbvXR92Pfv+pmn7v6WqdH1nTS371x48abzozdsGHDpP3gnbyuRwNpPhju6mvd3OHYrZ217iTVfLBbRpIqZLhLUoXsltGi0q3+bvvR1WsMdy0qnfR3dxLY9qOr1xju0hQMbPUz+9wlqUKGuyRVyHCXpAoZ7pJUIcNdkipkuEtShQx3SaqQ4S5JFTLcJalChrskVchwl6QKGe6SVCHDXZIq1LVwj4gdEXEmIoYjYn+31iNJullXwj0ilgB/Crwe2AK8JSK2dGNdkqSbdWvLfRswnJmPZ+YzwCFgZ5fWJUmaoFsX61gLPNl4PAq8ujlDROwF9paHT0fEmS7VMlPPB76z0EV0oN/qBWueD/1WL/RfzQtZ75TXcuxWuMckbTnuQeZB4GCX1j9rETGUmQMLXUe7+q1esOb50G/1Qv/V3Kv1dqtbZhRY33i8DrjQpXVJkiboVrh/DdgcEZsi4g5gN3CkS+uSJE3QlW6ZzLwWEb8FfB5YAnwsM091Y11d1LNdRlPot3rBmudDv9UL/VdzT9YbmTn9XJKkvuIZqpJUIcNdkmqUmdXcaB2h8yXgMeAU8I7SvgI4Cpwt98tL++3AIHCyLHOg8VqvKu3DwB9zowtrKfDJ0v4QsLGxzJ6yjrPAnlnWfE95/ENgYMIyB8r6zwDb57PmTusFfgk4Ueo6AbyuH37H5fmfAJ4GfqcfagZeCnylPH8SuLOH3xe9/Nn7MPAt4JvAA8DzeuGz1+ltwQN5Tn8YWAO8skw/F/h3WsMffAjYX9r3Ax8s028FDpXpHwFGxn75wHHgZ2kds/8PwOtL+28CHy3Tu4FPlukVwOPlfnmZXj6Lml8E/DTw5Qkfii3AN8qbZhPwbWDJfNU8g3pfAbygTL8Y+I/Gcz35O24s92ng7xgf7j1ZM62DI74JvKw8/vEef1/08mfvbuC20v5BbuTFgn72Or1V1S2TmRcz8+Ey/X1a38hraQ19MFhmGwR2jS0CLIuI24BnA88A34uINcCPZuZXsvWX+HhjmeZrfQq4KyIC2A4czcynMvM/af2HsGOmNWfmY5k52Vm7O2l9KK5m5jlaWwTb5qvmTuvNzK9n5tg5DqeAOyNiaY//jomIXbQ+cKcabb1c893ANzPzG2WZ72bm9V59X9Dbn70vZOa1MttXaZ2nM7b+BfvsdaqqcG+KiI20thofAlZn5kVo/UGBVWW2TwH/A1wEngD+IDOfovWFMNp4udHSBo2hFcob4L9pbSVNNuTCWjowoeapTLWeea+5zXqb3gx8PTOvLkS97dYcEcuAdwPvn/BUz9YM/BSQEfH5iHg4It61UDW3WW+/fPbeTmtLfNz6J6xnQWqeTreGH1hQEfEcWv9SvzMzv9f6opzUNuA68AJa/xr9S0R8kVsPnzDVc9MOudBJzbeadQbrn/OaO6h3bP6ttP7FvXuamm713Hz9jt8P/FFmPj3hvdPLNd8G/BzwM8APgGMRcQKYbJleeF/0/GcvIt4LXAM+MYv1d6XmdlS35R4Rt9P6Q30iMz9Tmi+Vf53G/rW+XNrfCnwuM/83My8D/wYM0PoWXdd42ebwCf8/tEL5l/LHgKeYxZALU9Q8lanWM281d1gvEbGO1o6pX8/Mbzdq6tXf8auBD0XECPBO4D3lpLxernkU+KfM/E5m/gD4e+CV81lzh/X29GcvIvYAbwJ+tXS1jFv/hPXMa81t67STvpdvtL4NPw58ZEL7hxm/Q/VDZfrdwF+W5ZYBp4GXlue+BryGGztI3lDa9zF+B8nhvLGD5BytrZDlZXrFTGtuPP9lxu+I2sr4nTqPc2OnTtdrnkG9zyv1vnmSeXvydzzhufcxfodqT9ZcXvthWjsnbwO+CLyxh98XPfvZo9X3fRpYOaF9QT97nd7mJXTn60br39KkddTAI+X2Blp9XMdoHXJ0bOyXCDyH1tEQp8of83cbrzUAPEprj/ifcOPQpjvLMsO09pD/ZGOZt5f2YeA3Zlnzr9D6dr8KXAI+31jmvaWuM5S98vNVc6f1Ar9Hq2/1kcZtVa//jhvLvo/x4d6zNQNvo/VefpSyAdPD74te/uwN0+oPH2v7aC989jq9OfyAJFWouj53SZLhLklVMtwlqUKGuyRVyHCXpAoZ7pJUIcNdkir0fz4pepUS5JjJAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Zoom in near the miniumum to see the quadratic shape near the best-fit value w0Star.\n", "plt.plot(w0List, ChiSqList, 'ks', linewidth = 2.5, markersize = 5, fillstyle = 'none')\n", "plt.axis((w0Star - 15*w0step, w0Star + 15*w0step, 0, 500));" ] }, { "cell_type": "code", "execution_count": 68, "id": "floral-pioneer", "metadata": {}, "outputs": [], "source": [ "# The last step is to characterize the shape of the quadratic dependence:\n", "\n", "# ChiSq = ChiSqZero + B*(w0 - w0Star)**2\n", "\n", "# There are three unknowns (ChiSqZero, B, and w0Star). These could be determined\n", "# by fitting the data above to a quadratic finction. Alternatively, if we pick\n", "# three values of w0 near the minimum, we can use the corresponding ChiSq values\n", "# determine the unknowns.\n", "ChiSq1 = ChiSqList[244]\n", "ChiSq2 = ChiSqList[247]\n", "ChiSq3 = ChiSqList[250]\n", "dw0 = w0List[250] - w0List[247]" ] }, { "cell_type": "code", "execution_count": 69, "id": "characteristic-booth", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "w0min = 213865.68039978252\n", "Δw0 = 304.59089220406986\n", "ChiSqZero = 48.14426447267113\n" ] } ], "source": [ "# The value of w0 that minimizes ChiSq (w0min) is given below as is the\n", "# uncertainty in w0min which is determined by how quickly ChiSq increases\n", "# as we move away from w0min. That is, errw0 is determined by the B\n", "# coefficient in the quadratic depedence of ChiSq on w0. The larger the\n", "# value of B, the lower the errw0.\n", "w0min = w0List[250] - dw0*((ChiSq3 - ChiSq2)/(ChiSq1 - 2*ChiSq2 + ChiSq3) + 1/2)\n", "errw0 = np.sqrt(2)*dw0/np.sqrt(ChiSq1 - 2*ChiSq2 + ChiSq3)\n", "ChiSqZeroFit = ChiSq3 - (w0List[250] - w0min)**2/errw0**2\n", "print('w0min =', w0min)\n", "print('\\u0394w0 =', errw0)\n", "print('ChiSqZero =', ChiSqZeroFit)" ] }, { "cell_type": "code", "execution_count": null, "id": "endangered-shame", "metadata": {}, "outputs": [], "source": [ "# Notice that the extracted value for the uncertainty in w0 is very similar to the value \n", "# reported by the 'curve_fit' routine from the scipy.optimize module.\n", "\n", "# In general, extracting uncertainties in the fit parameters from a nonlinwar fit is\n", "# nontrivial. Above we have given one method for obtaining an estimate if the \n", "# uncertainties in the fit parameters. For more details about uncertainties for\n", "# nonlinear fits, see Data Reduction and Error Analysis by Bevington and Robinson\n", "# and/or Numerical Receipes in C which you cna obtain and read fro free by visiting\n", "# http://www.nrbook.com/a/bookcpdf.php." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.8" } }, "nbformat": 4, "nbformat_minor": 5 }